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INTRODUCTION 


Work supported under Contract Number NAS8-28248 began In 
January 1972 and has continued with varying degrees of effort to the 
present time. The work has Involved the application of improved aero- 
dynamic theory towards the goal of obtaining more accurate knowledge 
of the upper atmosphere. Advancements have been made both in the area 
of aerodynamic theory and the interpretation of the dynamic response 
of objects traveling through the atmosphere. 

The work began as a study of ways of improving models of the 
upper atmosphere as deduced from observation of satellite decay. In 
the development of atmospheric models, it was found that the decay 
of the orbit of a satellite due to drag had been modeled as simply a 
sphere with a drag coefficient of 2.2 traveling through a rotating 
atmosphere. The assumption of a sphere of C^ = 2.2 was of course 
recognized to be only an approximation and of particular use in the 
analysis of the drag decay of foreign or some domestic satellites for 
which limited knowledge existed concerning the shape or other physical 
parameters. One of the goals of this work was to investigate the 
magnitude of error made in the assumption of a sphere of C = 2.2 and 
to propose more accurate data reduction techniques. 

Chapter 1 describes the major influence revealed in this study 
concerning the influence of real satellite aerodynamics on the deter- 
mination of upper atmospheric density. Chapter 2 presents a method 
of analysis of satellite drag data which includes the effect of satellite 
lift and the variation in aerodynamic properties around the orbit. 


The method was applied to the data from OVI 15 satellite. One of the 
interesting results obtained from analysis of satellite orbit decay 
has been the super rotation of the atmosphere deduced by King-Hele. 

In Chapter 3, a study is presented which shows that satellite lift 
effects may be responsible for the observed orbit precession rather 
than a super rotation of the upper atmosphere. 

Emphasis of this work gradually came to the lower altitude regime 
and the 80 to 120 km region in particular. This region of the atmo- 
sphere is of great Importance since it serves as a major boundary 
between the lower atmosphere which is constant in molecular weight and 
the upper atmosphere which reaches out many thousands of kaj having 
considerable variation in molecular composition. This Important region 
of the atmosphere has received little experimental attention due to 

X i ' ' 

the difficulty and expense of performing measurements at these altitudes. 
The falling sphere method is found to be the primary source of information 
for this region. As with satellite drag measurements, it was found 
that simple assumptions concerning the aerodynamics of objects were 
often employed in the falling sphere analysis. The Influence of the 
errors made due to the simplifying assumptions were evaluated and an 
Improved method of analysis was proposed and applied as reported in 
Chapter 4. 

The work on falling sphere data analysis also revealed that 
most of the 80-120 k.m results were based on values of drag coefficient 
in the transition regime that were extrapolated from wind tunnel results 
that were far outside the transition regime. In the work reported here, 
more recent wind tunnel data reported in the literature were obtained 
and more accurate drag coefficient relationships were developed based 
on these data. This work is reported in Chapter 5. 
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The improved drag coefficient relationships revealed a con- 
siderable error in previous falling sphere drag interpretation. These 
data were reanalyzed using the more accurate relationships. The 
result of this work is given in Chapter 6. 

In this work the drag coefficient has been studied for the 
entire spectrum of Knudsen Number and speed ratio. One region which 
was of particular interest is in the very low speed ratio region. 

This region of the aerodynamic spectrum has received little experimental 
attention except for highly viscous flows due to the experimental 
difficulty of obtaining drag data jin a low-densitjj slow flow situation. 
The theoretical work in this region is discussed in Chapter 7. 

The recommendations for future work are presented in the form 
of two proposals given in Chapter 8. Both proposals would involve 
additional analytic work and subsequent experiments using the shuttle 
space craft. 

The computer programs generated during the period of performance 
of the contract are given In the appendix. 


CHAPTER I 


INFLUENCE OF SATELLITE AERODYNAMICS ON ATMOSPHERIC DENSITY DETERMINATION 

A preprint of a paper by G. R. Karr and R. E. Smith from Preprint 
Volume of the International Conference on Aerospace and Aeronautical 
Meteorology was published by AMS, Boston, Mass. Presentation was given 
May 22-26, 1972. 
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RaprlaUd (rom Pniirliit Vohimt of tha htonatlonal Canforoaea on Aanapaca and Aaronaatleal 
Mataorototr, May 11-16, ISTl, Waahli«taa, D.C.; publlahad by AW, Boaton, Maas. 


INFLUENCI or SATELLITE AEttOOnUiaCS ON ATMOSnEUC DENSITY DETERMINATION'*' 


Gerald R. Karr 

Coordinated Science Laboratory 
Unlveralty of Illlnola 
Urbana, Illlnola 


1. INTRODUCTION 

Drag-deduced denaltlaa of the upper atmoa- 
phere have been a primary aourca of data In the 
development of atmoapharlc modela and to the 
atudy of the upper atmoaphere. In the paat, the 
determination of atmoapharlc denalty haa been 
through the obaarvatlon of aatelllte orbital 
decay ovar a long period of time which necea- 
aarlly required knowledge of only the average 
drag propertlea of tha aatalllte. Howevar, aa 
Cracking technlquea become more accurate and the 
uae of aenalclva accelerometera Increaaaa, the 
aaaumpclon of average drag propartlaa la no 
longer valid and a more accurate treatoMnC of 
aatelllte aerodynamics oaiat be awda. 

The purpose of Che following discussion will 
be to focus on three principle saCelllta aero- 
dynamic factors which Influence tha Interpre- 
tation of aatelllte dynamic response; these are, 
(1) the Influence of satellite orientation and 
shape on the drag coefficient, (2) the effect of 
changes In tha gas flow propertlas with altitude, 
and (3) the Influence of upper atmospheric winds 
on the IntarprsCaClon of data. 

The three topics to be treated are effects 
causing tha greatest source of error In current 
data reduction. Othar factors such as aero- 
dynamic lift, changing atmospheric composition, 
and changing aatalllte surface propertlea will 
not bo traated hara but such factors could be of 
Importance for particular satellite systems 
having large aerodynamic lift forces and wldaly 
varying gas and surface propertlas. The follow- 
ing will than ba limited to a discussion of 
aerodynaadc drag affects only and the assumption 
of constant satelllta surfaco properties. Tha 
atmospheric gas will be considered of single 
species having the average properties assoclatad 
with a particular altitude. 


This work was supported In part by tha Joint 
Services Electronics Program (U.S. Army, U.S. 
Navy, and U.S. Air Forca) under Contract DAAB- 
07-67-C-0199; and In part by tha National 
Aeronautics and Space Administration and the 
American Society for Englnsaring Education 
through the 1971 ASEE-NASA Sunmiar Faculty 
Fellowship Program at Marshall Space Flight 
Center . 


and Robert E. Smith 

National Aeronautics and Space Adadnlatratlon 
Marshall Space Flight Center 
Huntsville, Alabama 


The three factors to be discussed are of 
Importance In the Interpretation of past data as 
wall as the application to future more accurate 
measureawnts . For this reason, an estimate will 
be siada of the possible correction to present 
density mdels based on the results of the 
present study. At a basis for calculation and 
comparison only, the Jacchla (1971) s»del Is 
employed in the analysis. Other current models 
could be used for the same purposes with similar 
results . 

2. SATELLITE AERODYNAMICS 

The aarodynaaiic flow regime experienced by 
the majority of satellites Is free molecular. 

That la, tha mean free path between collisions 
of molaculat In the upper atmosphere Is greater 
then the dimensions of the majority of Earth 
satellltaa. This assumption Is certainly true 
for altitudes greater than 200 km where the near 
free path Is on the order of kilometers. 

Although departures from free molecular flow will 
occur In regions of high density, the assumption 
of free s»lecular flow may be considered 
reasonably accurate to an altitude of 100 km 
where the mean free path Is of the order of 
meters. For convenience, free molecular flow 
la assumed through out the altitude range being 
considered in this discussion. 

In the free moleculer flow of space, by 
definition, collisions of molecules with the 
satellite surface predominate. For this reason, 
the atudy of satellite aerodynemlcs requires an 
understanding of tha Interaction of gas sxilecules 
with solid surfaces. The drag properties of a 
aatalllte are Influancad prlsiarlly by tha 
axchangs of momentum with the surface during the 
Bulecular Impact. 

The description of the swilecular Impact to be 
employed in this discussion Is the generalized 
gas surface interaction (GSI) modal which has 
direct application to tha problem of satellite 
dreg, sea Karr (1969). This modal Is represented 
in Figure 1 shoi^ng a general non-specular type 
of reflection. The reflected molecules produce a 
mofsantum vector In the direction 8. with an 
average velocity of U, . Tha reflected properties 
are assumed to depend^upon the Incident flow 
properties such that 
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9j - 2 Pj + a-fj)9 


Cjj - 2 -zyi^ ■*• (2 -Pj)6] 



Fig. 1. Description of the gas surface inter- 
action . 


Using the model of the interaction described, 
the force components In the direction of drag can 
be expressed locally at the surface. The total 
force is obtained by integrating over the entire 
surface exposed to the flow. If the satellite 
velocity is much higher than the random kinetic 
motion of the gas molecules, the gas molecules 
can be considered stationary as the satellite 
sweeps out molecules in its path. The assumption 
of such conditions (the hypervelocity assumption) 
allows one to determine the drag coefficient 
based upon the area projected to the flow. 

Under the assumption of hyperveloclty flow 
and Che generalized gas surface interaction Che 
following results are obtained for four shapes of 
Interest. The drag coefficient is defined as 

Cjj ■ Drag/^P V^A 

where A is a reference area taken to be a 
constant, independent of angle of attack. 

Flat plate with angle of attack 9 

A “ Area of plate 

Cjj “ 2 slnS - 2Vl-»j sin9 cos J Pj + (2-Pj)P 


Cylinder with axis perpendicular to flow 



length 


where the bracketed term is equal to ^ when P^ ” 

1,0. A value of P - 3 is not meaningfully applied 
to this shape. ^ 


Cone with axis parallel to flow 
A " Area of base of cone « n r^ 


where 6 is Che cone half angle. 


Sphere 


A 


n 


2 

r 



2 + 



■ TT * 

4(l-cos - Pj) 

P (4-P ) 


where the bracketed term is equal to zero for 
Pj ■ 0. A value of P^ “ 4 is not applicable. 

The above results clearly Illustrate Che 
influence of Che GSI and the shape on the drag 
properties. The results for a sphere are 
shown plotted in Figure 2 . A range of Cjj from a 
minloaim of 2.0 Co a maximum of 4.0 is seen 
depending upon the values of gas surface inter- 
action parameters, or^ and . The results for 

the cone shape reveals that depending upon Che 
cone half angle, Cq values less Chan one are 
possible. Thus, it is seen that satellite 
shape and satellite surface properties are 
strong Influences on the drag properties. 



Fig. 2. Drag coefficient of sphere as a function 
of gas surface interaction parameters. 


3. THE EFFECT OF SATELLITE ORIENTATION 

Satellites rarely present the same shape to 
the flow during an orbit. For this reason, the 
effect of satellite orientation, which is 
important in determining the instantaneous drag 
properties of a satellite, is also important in 
determining the average drag properties. In 
order to illustrate the effect of angle of 
attack, consider first a cylinder with spherical 
end caps with an angle of attack Ot. For the 
special case of specular reflection (pj>0) and 
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hypervelocity flow, the following analytic 
results are obtained. 


Cylinder with spherical ends ; Pj =0, hyper- 
velocity flow and angle of attack Qf. 

Cd = 2 [l 1^5 sino] 

+ '/l-a . ^ i sin siii^CT - 2.0^ 

_ 2 

where A = and l./D in tlie length to diameter 
ratio of the cylinder. 

These results reveal that the cylinder 

with spherical ends, changes considerably depend- 
ing upon the angle of attack and the shape 
(length to diameter ratio). As expected, the Cjj 
value is that of a sphere when o is zero and is 
progressively influenced by the cylinder as the 
angle of attack is increased to a value of 

j. Since n r^ is used as the reference area at 

all angles of attack, will reach high values 
for large values of L/D. 


For more complex shapes such as the cone, 
results must bo obtained numerically. Figure 3 
shows such results for a 35° half angle cone. 

The plot is a three axis presentation of the 
surface Cq(P.,P) where P. varies from zero at 
the front to^a value of at the rear. The 
angle of attack is varied from zero on the right 
to 180° on the left. Tlie value of a. was taken 
to be zero in constructing this plot. The 
function has a maximum value of 4.00 and a 
minimum value of 0.399. 'Ihese results illustrate 
further the strong dependence of aerodynamic 
drag on the shape, orientation, and CSI . 



o 


For a circular orbit the satellite will travel at 
a constant rate and, for this case, 



Consider now a spin - stabilized cone shaped 
satellite with a flat base. The orientation of 
the satellite spin axis with respect to the 
orbit Is given by the angle X as shown In 
Figure 4. The angle 0^ serves to orlente the 




Fig. 3. Computer generated plot of surface 

Cq(P., 3) for a 35° half angle cone with 

flat base where 1'. varies from zero to 
two from front to^back and 0 varies 
from zero to 180” from right to left. 
The values of vary from 0.399 to 
4.00. 

Just as the instantaneous values of Cq are 
Influenced by angle of attack, the average value 
of Cjj over an orbit also depends upon the 
attitude history the satellite experiences over 
the orbit. in order to illustrate this factor, 
consider the determination of the average drag 
coefficient over an orbit of the satellite. The 
time average of C^ is given by 


Fig. 4. Coordinated system describing 

orientation of spin stabilized satellite 
with respect to orbit plane. 


£pln axis with respect to the velocity vector 
U . For a circular orbit, the velocity vector 
will rotate about the spin stabilized satellite 
at a constant rate. Since the instantaneous 
value of Cq for each point in the orbit can be 
determined, the average Cq over the orbit can 
also be determined by numerical quadrature. The 
results of such a study are presented In Figure 
S for five cone half angles and for X from 0 to 
90° . The base area of the cone was chosen as the 
reference area in each case and the same GSl 
parameters used for each plot. 

The results given In Figure S Illustrate the 
Importance of taking the satellite orientation 
history into account In selecting an average 
drag coefficient. Since the satellite orbit and 
the satellite spin axis will tend to drift In 
space, the average drag coefficient should not be 
expected to remain constant. The amount of 
variation Is seen to Increase as the amount of 
non-symmetry of the satellite Is Increased. Only, 
spherically synmetric shapes will experience no 
variation. 
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Fig. 5. Average drag coefficient for apln 

stabilized cone shaped satellites as a 
function of spin axis orientation with 
respect to orbit. 

4. EFFECT OF ALTITUDE ON DRAG COEFFICIENT (THE 
SPEED RATIO EFFECT) 

The assumption of hypervelocity flow used in 
obtaining the results of the preceding sections 
will now be examined. A discussion of this 
assumption Is facilitated by Introducing a flow 
parameter termed the speed ratio, S, defined as 
the ratio of the satellite velocity to the 
random velocity of the gas molecules In thermal 
equilibrium. The speed ratio Is given by 



Fig. 6. Circular speed ratio as a function of 
altitude for four exospheric tempera- 
tures based on Jacchla 1971 model. 


case. These three values were used to obtain a 



S « U^/72RT/M 

where R the universal gas constant, M the 
.IV rage molecular weight and T Is the tempera- 
ture of the gas. Hypervelocity flow, the 
assumed flow In the preceding work. Is approached 
as the speed ratio approaches Infinity. 

Two factors Influence the value of the speed 
ratio as the altitude of the satellite orbit Is 
changed. First, considering circular orbits 
only, the velocity of the satellite decreases as 
altitude Increases. Second, the temperature of 
the atmosphere and therefore the thermal velocity 
of the molecules Increases with Increasing 
altitude. These two factors combine to cause a 
considerable decrease In speed ratio as altitude 
Increases. This effect Is Illustrated In 
Figure 6 which shows the value of the circular 
speed ratio aa a function of altitude for four 
exospheric temperstures . This plot was 
constructed using values of temperature and 
mean molecular weight from Jacchla (1971). 


• 2.0 

8 2.106 + 0.450(2L/ttd) 

4 2.249 + 0.885(2L/ttd) 


second order polynomial approximation to the 
variation In Cg with S. The results are 

C /C - 1 + (0.350 + 1.166 L/D)/S 

S- 2 

+ (0.592 - .1528 L/D)/S . 

This function Is plotted In Figure 7 for three 
values of L/D. Two factors of Importance are 
Illustrated In these results. First . long 
slender object at low angles of attack are 
strongly Influenced by speed ratio effects. 

Even at relatively high values of speed ratio, 
there Is strong sensitivity to the length to 
diameter ratio for long slender objects. 

Second , for a given satellite shape, consider- 
able change In the value of Cq Is seen to occur 
over the range of Interest from S°4 to S=20. 


In order to Investigate the Influence of 
speed ratio on the drag coefficient, consider a 
cylinder with spherical ends. Utilizing data 
from Karr and Yen (1972), Sentman (1961), and 
Fan and Andrews (1969), the value of for 
three speed ratios were obtained giving the 
following results for the zero angle of attack 


Combining the results of the polynomial fit 
of Cq values for a sphere with the change of S 
with respect to altitude, the change In Cq for a 
sphere over the altitude range from 100 to 1000 
km Is obtained and presented In Figure 8. These 
results, for the same four exospheric 
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S. EFFECT OF ATMOSPHERIC HINDS ON DRAG STUDIES 



Fig. 7. Change in drag coefficient as a function 
of speed ratio for three values of L/D. 



Fig. 8. Change in drag coefficient as a function 
of altitude for four exospheric 
temperatures . 


temperatures as In Figure 6, show that a sphere 
will experience a change In Cq of from 2 to 107. 
depending upon the atmospheric temperature. 

Even greater changes would be expected for non* 
spherical satellites such as the cylinder with 
spherical ends. 


Recent analytic and experlsMntal studies 
point to high velocity atsnspherlc winds In the 
upper atsK>sphere. Velocities of as nuch as 
400 m/sac have been reported. Since satellite 
velocities are of the order 7 lus/sec, atmospheric 
winds are expected to be an Influence on the 
satellite ax>tlon. In order to Investigate such 
Influence, consider a simple model of the upper 
atsK>spheric wind structure which Includes the 
rotation of the atmosphere at the rotation rate 
of the Earth. In addition to the rotation, 
consider an east-west wind which varies In both 
magnitude and direction. Results of Challlnor 
(1969), suggest that the east-west component 
could be assumed to approximate a sinusoid 
variation with a longitude angle, or, measured 
from a reference .point of zero wind. The 
velocity of a satellite with respect to the 
atmospheric gas in circular orbit Is then given 
by 

U “ tCi -V sina 

to r e max 

where V Is the peak wind velocity, r Is the 
dlstance*^from the center of the earth and 0^ 

Is the angular velocity of the earth. At some 
longitudes, the wind Is seen to subtract from 
the atmospheric velocity while adding at other 
longitudes . 

The Instantaneous drag acting on a satellite 
at any point In the orbit Is given by 

D - ^ PU^CjjA. 

The average drag over one revolution is found by 
Integrating over the time In Che circular orbit 

- J PC A I j 1 

the average drag Is found to be increased by a 
factor, where 

^ max 

•^wlnd ■ 

r e 

The Increase In drag due to winds for a circular 
orbit Is found to be small enough to be 
neglected. The Increase Is less than 17. for 
typical values of velocities. 


Of importance is the fact that the noted 
change in C- with altitude Is found to be 
systematic with altitude. For this reason then, 
it Is to be expected that present drag deduced 
density values will have these systematic errors 
incorporated Into the results. This factor 
could account for some of the discrepancies in 
densities at high altitude in comparison to 
densities at lower altitudes. Similar results 
and conclusions were found by Izakov (1965) 
using an analytic expression for Cq of a sphere 
as a function of S . 


Consider now a more severe case when a 
satellite In an elliptic orbit has Its perigee 
at the peaks of wind velocity. Such a situation 
Is likely since, often tlsws, observations of 
elliptic satellite orbits are made to determine 
atmospheric properties In the region of perigee. 
Much of the knowledge of latitude and longitude 
variations In the atmosphere have developed 
from such observations. 

Consider the perigee passage at ar» 90°, 
where the wind velocity subtracts, with the 
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perigee peeeage at 270**, where the wind velocity 
adda. For aiailar orbita and conatant denaity, 
the ratio of drag at perigee for these two eases 
could be taken to 


"90° 
°270° ■ 


./^^/^+e - r n 


P a 


yi +e - r n + 


V 

p e MX 


4 V 




vnr - r n 

r p e 


where r is the orbit radius at perigee and e is 
the eccentricity of the orbit. In the altitude 
range from 100 to 1000 km, r^C)^ is approxiiaately 

.5 km/sec, and Vti/rp is approximately 7.5 km/aec. 
For e values less than 0.2, percentage 
differences in drag of Che order of three per- 
cent are obtained for values of 200 to 

400 m/sec. The wind effect is found to be such 
larger for this case than for Che circular orbit 
case. Due to the lag in solar heating of the 
atmosphere, the 0-90° point would occur after 
sunset while Che 270° point would be after sun- 
rise. The results obtained reveal that a 
difference in drag of a maximum of 3% at about 
200 km is explainable by wind effects. If 
these wind effects were not taken into considera- 
tion, the difference in drag would be misinter- 
preted as being caused by a corresponding 
difference in atmospheric density. Therefore, 
wind effects could explain some of the day-night 
variation in density deduced from satellite 
drag. 


6. CONCLUSIONS AND DISCUSSION OF RESULTS 

The preceding discussion has emphasized the 
variability of satellite drag coef ficianCa . In 
particular, Che effect of the gas surface 
interaction is seen to dominate. Unfortunately, 
values of P. and or have yet to be determined 
accurately enough to be used in satellite drag 
studies. This factor leads to considerable 
uncertainty in the specification of satellite 
drag coefficients which could cause errors of 
as much of 50% in the values of Cq ■* 2.2 used 
in past data reductions. 

In view of the influence on drag coefficients 
of non-spherica lly symnetrlc shapes and the 
speed ratio, values of Cq of around 2.2 are 
likely Coo low for most satellite shapes. This 
observation is based on the fact that for a 
sphere Che minimum Cq for S-** is 2.0. Speed 
ratio effects cause the minimum value to be 
increased to values of 2.1 or higher. In order 
for the sphere Co have the minimum Cq, the gas 
surface interaction would have to be specular 
or the accommodation coefficient would have Co 
be unity correaponding Co reflection at near 
zero velocity. Such limiting interactions 
appear unlikely, meaning that Cq values highar 
than 2.2 are likely. 

Although an absolute value for Cq may be 
lacking presently, the results obtained for 


changes of Cq dua to shape, orientation, and 
altitude reveal that systematic variation of Cq 
around an absolute value would occur under most 
circumstances. These variations should be 
Included In the analysis of future satellite 
data where Inforsmtlon on satellite orientation 
may be available. 

The effects of atmospheric winds were 
Illustrated assuming a simple model of the wind 
structure. These results Illustrate still 
mtother source of error in drag deduced 
density values. Future planned work will 
consider a more sophisticated model of the wind 
structure. This study Is expected to lead to 
Improved knowledge of the upper atmospheric 
density and wind structure. 

7 . REFERENCES 

Challlnor, R. A., 1969, Neutral-air Winds in the 
Ionospheric F-raglon for an Asymmetric 
Global Pressure System, J . Atmos . Scl . . 17, 
pp. 1097-1106. 

Fan, Chlen and Andrews, C. Donals, 1969, 

Application of the Nocllla Wall Reflection 
Model to Calculate Aerodynasde Coefficients 
for Orbital Vehicles with Complex Geometry, 
Lockheed Missiles and Space Company, 
HBEC-4401-1, Huntsville, Alabama. 

Izakov,' M. N., 1965, Some Problems of 

Investigating the Structure of the Upper 
Atmosphere and Constructing its Model, Space 
Res., V, pp. 1191-1213. 

Jacchla, L G., 1971, Revised Static Models of 
Che Thermosphere and Exosphere with 
Empirical Temperature Profiles, SAO Special 
Report 332, Smithsonian Institution, 
Casdirldge, Mass. 

Karr, Gerald R., 1969, A Study of Effects of the 
Gas-Surface Interaction on Spinning Convex 
Bodies with Application to Satellite 
Experiments, Coordinated Science Laboratory 
Report R-43S, University of Illinois, Urbana, 
Illinois (Rt.D. thesis). 

Karr, G. R. and Yen, S-M, 1972, Aerodynamic 
Properties of Spinning Convex Bodies In a 
Free Molecule Flow, Rarefied Gas Dynamics . 
Seventh Symposium, Academic Press, (In 
press), 8 pages. 

Sentman, Lae H., 1969, Free Molecule Flow Theory 
and Its Application to Che Determination of 
Aerodynamic Forces, Lockheed Missiles and 
Space Tech. Report, LMSC-448514, Sunnyvale, 
California. 


10 



:b 


CHAPTER II 

SATELLITE AERODYNAMICS AND ATMOSPHERIC DENSITY DETERMINATION FROM 

SATELLITE DYNAMIC RESPONSE 


11 



SATELLITE AERODYNAMICS AND ATMOSPHERIC DENSITY DETERiMlN^TI 0 N FROM 

SATELLITE DYNAMIC RESPONSE 

by 

Gerald R. Karr 


ABSTRACT 


A method for determining satellite aerodynamic properties and upper 
atmospheric density from obsc-rved satellite dynamic response has been 
successfully developed and tested. 

The aerodynamic drag and lift properties of a satellite are first 
expressed as a function of two parameters associated with gas-surface 
interaction at the satellite surface. The dynamic response of the 
satellite as it passes through the atmosphere is then expressed as 
a function of the tv/o gas-surface interaction parameters, the 
atmospheric density, the satellite velocity, and the satellite 
orientation to the high speed flow. By proper correlation of the 
observed dynamic response with the changing angle of attack of the 
satellite, it is found that the two unknown gas-surface interaction 
parameters can be determined. Once the gas-surface interaction 
parameters arc known, the aerodynamic properties of the satellite at 
all angles of attack are also determined. The atmospheric density 
may then be accurately calculated once the true aerodynamic properties 
are knovm. 

Employing accelerometer data from the OVl-15 satellite, analysis was 
successfully made of the aerodynamic properties of that satellite 
and a determination was made of the absolute value of atmospheric 
density near the orbit perigee. These results constitute the first 
successful application of the proposed method of analysis. These 
results also serve to illustrate the potential of the technique in 
the analysis and prediction of satellite orbit decay in the atmos- 
phere and the accurate determination of upper atmospheric density 
from satellite dynamic response. 
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1 . Introduction 


The problem of satellite orbit decay prediction and the problem of 
upper atmospheric density determination have encountered a common 
source of unknown error which can be traced to a lack of knowledge of 
satellite aerodynamics. The basic equation employed In both orbit 
decay and density determination is the familiar drag equation: 

Drag = i/2 pU^Cp)A 

where p is the density, U is the velocity of the satellite with respect 
to the atmosphere, Cj) is the drag coefficient, and A is a suitable 
reference area. In most applications of this equation to sateHltes 
the value of Cj^ is considered to have a constant value. Generally, 
however, the assumption of constant drag coefficient is not valid and 
the use of such an assumption can lead to considerable error (see 
Karr, 1972). The uncertainty in satellite aerodynamics has prevented 
the assignment of even an approximate value of drag coefficient with a 
known range of uncertainty. A value of Cd of 2. 0 or 2.2 is often used 
in satellite drag studies and in the determination of atmospheric density. 
These values of Cp) are likely too small and, combined with the fact 
that Cj) is not constant, have resulted in an overestimation of upper 
atmospheric densities (see Karr and Smith 1972), 

A more accurate treatment of satellite aerodynamics has obvious 
benefit to the determination of upper atmospheric density and the pre- 
diction of satellite orbit decay. Satellites traveling in the earths upper 
atmosphere experience the aerodynamic flow regime termed the free 
molecular flow regime. In this flow regime, the collision of atmospheric 
gas molecules with the satellite surface dominate the flow and collisions 
of gas molecules with other gas molecules may be ignored. Satellite 
aerodynamic properties are then the result of the interaction of high 
speed gas molecules with the solid satellite surface. Unfortunately, 
very little is known about the gas surface interaction at satellite veloci- 
ties and this lack of information is the basic source of uncertainty in 
satellite aerodynamic properties. 

In the interest of developing a more accurate treatment of satellite 
aerodynamics for application to orbit decay prediction and density determ 
nation, a model of the gas surface interaction has been developed which 
utilizes two parameters to describe the interaction (sec Karr, 1969 and 
Karr and Yen, 1970). The advantage in this treatment of satellite 
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aerodynamics is that no a priori assumptions of the aerodynamic / 
properties need be made. The gas surface interaction parameters 
are considered as unknowns to be determined from the observed / 
dynamic response of the satellite as it travels through the atmosphere. 
The determination of the gas surface interaction parameters serves 
as the key to the subsequent determination of both the aerodynamic 
properties and the atmospheric density. 

To test the proposed method of analysis, accelerometer data from 
the OVT-15 satellite is used. The accelerometer data provide an 
accurate, instantaneous measure of the level of aerodynamic force 
and the attitude of the satellite with respect to the flow. As is pointed 
out in the paper, accurate satellite attitude informatin is essential to 
the analysis. The OVI-15 satellite, although less than ideal in shape 
for an aerodynamic study, provided a good basis for the test of the 
proposed method of analysis. The results serve to illustrate the 
potential that this method of analysis has to future determinations of 
aerodynamic and atmospheric properties. 

2. Satellite aerodynamics and the gas surface interSi.ction. 

Consider a local satellite surface element in which the high speed 
flow of molecules is incident at an angle of 0 as shown in Figure 1, 
As.socLated with the incident flow is the incident moinentum which 
gives rise to the incident force This force is colinear with the 

satellite velocity, U, with respect to the atmosphere. Assume for 
now that the speed ratio is infinite where the speed ratio is defined 
as the satellite velocity divided by the thermal velocity of the gas 
molecules. The thermal velocity of the gas molecules is taken to be 
equal to 4 R T/M where R is the gas constant, T is the temperature, 

M is tJie mean molecular weight. 

The molecules reflected from the surface cause a net reaction 
force Fj. which is colinear witJi the mass-motion velocity vector Uj 
of the molecules leaving the surface. The direction of Uj is given 
by the angle 0^. 

Modeling of the interaction is performed by providing relationships 
between the incident and reflected quantities. The relationships are 
given by 

U- 1 - a . U 

J J 

0j= ^ Pj + (1-pj) e 
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where Of j and Pj are the parameters of the interaction. The subscript j 
is used if more than one set of interaction coefficients is to be considered 
in the analysis. Vnbirl more is learned of the interaction, these linear 
relationships provide a useful first approximation to the interaction that 
occurs at satellite velocities. The parameters Ofj and Pj ax’e capable 
of describing a much wider range of possible interactions than other 
models. The development of this model and the capabilities are 
described in detail in Karr 1969 and Karr and Yen, 1970. 

The model described above is particularly useful in the determini- 
nation of forces acting on the satellite surface. The total vector force 
acting on the element of surface shown in Figure 1 is given by 

dF = -(U-CTjUj) p U-n dA 

where CTj is employed if more than one gas surface interaction is 
employed in the analysis. In order to conserve mass at the surface, 
the sum of the a j values must be unity. 

The magnitude and direction of Uj is determined by the parameters 
Ofj and Pj. The vector force acting on the local element of surface is 
then expressed as a function of OCj, Pj, U, p, and 0 . For a given satellite 
shape (assumed to be convex^ the total aerodynamic forces and torque 
acting on the satellite are found by integration of dF and RxdF over the 
surface exposed to the flow. In general, the results will be of the form 

Drag = 1/2 p Cj 3 (Ctj. Pj, P) A 

Lifti 2 = 1/2 p U2 ^ 

Torque^ 2.3 = l/ZPU^ ^ Pj. P) AX 

where 0 is an angle of orientation and the subscripts on Cj^ and C-p are 
to indicate that there are two coinponents of lift and three components 
of torque. The six aerodynamic properties are found to be a strong 
function of the gas surface interaction parameters. For non-spherical 
objects, the angle of orientation, P , also has a strong influence on the 
drag, lift and torque properties (sec Karr and Yen, 1970). 
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3. Aerodynamics of the OVI-15 satellite. 

The approximated shape of the OVT-15 satellite is a cylinder with 
spherical ends as shown in Figure 2. The satellite spin axis was 
normal to the longitudinal axis of the cylinder with a spin rate of about 
10 rpm. Near the center of the satellite a three axis accelerometer 
detected the forces acting on the satellite. Since the data to be used in 
the subsequent analysis has been filtered and averaged over a number 
of spin cycles,. the aerodynamic properties averaged over a spin cycle 
are developed. 


The total instantaneous vector force acting on the satellite is given 
by 


p = 1/2 P K {Gjj D F Clj Li + Cl2 Lz) 


where Dj Ljj and L£ mutually orthogonal unit vectors in the drag 
and lift directions. The Lx and L 2 directions are defined with respect 
to the instantaneous orientation such that Lj is perpendicular to the 
cylinder axis and the velocity vector. The direction of L 2 perpendi- 
cular to both Lj and D. Due to symmetry the lift force in the Li 
direction is zero. 


From Karr, 1969, Cq and Cj ^2 obtained for the infinite speed 
ratio case, given by 


Cj 3 = 2 + 4 Vl -«j (1 -cos ^1^) Pj (4 - Pj) 


+ 2 Ap^ cos 0, 


2TT 


+ Ap V 1 -0^ Cj - cos3 


sin^ § (Cj !-Sj) ]d 5 


2it 





where the first two terms in Cj) are due to the spherical ends and the 
remaining terms are due to the cylindrical section. The angle ^ is a 
cylindrical surface-integration angle. The quantity Ap is the area ratio 
of the cylinder to the sphere given by 

Ap = ^cyl/-^sph “ 2rL/^^r = 4 l/tt D 
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The quantities C; and S- contain the parameter P- where 
J J J 

COS (l r/2 Pj -I- (1 - Pi) 9) 


Cj = 


Sj = 


cos 9 

sine/2 Pj +(1-Pj) 6) 


sin 0 


-1 


0 = sin” (— cos 0g sin 5) 


The angle 0^ is an angle of instantaneous orientation defined as the angle 
between the velocity vector and the longitudinal axis of the cylinder. 

Since the OVI-15 spin axis is perpendicular to the axis of the cylinder, 
the angle 0 s is a function of the angle the spin axis makes with the velocity 
vector, Y» and a spin angle, , which changes from 0 to 2n every spin 
cycle (see Figure 3). 

For certain values of Pj the surface integrals over the angle 5 are 
easily performed. For Pj = 0, which corresponds to specular type 
reflection. 


= S 


P-O 


^|pj=0 


= 1 


For Pj = 1 which corresponds to diffusive type reflection, 




= 0 


Pj=i 




Pj=i 


= -l/cos 0 sin § 


For Pj = 2 which corresponds to perfect backscatter type reflections 


= -1 


Pj=2 


= 1 


P.= 2 


The values of Ct^ and Cj_,2 three values of P- = 0, 1, 2, were 

used to obtain an polynomial approximation for Cp> and Ct as a function 

of Pj. 

since the accelerometers were body fixed, the outpiit of the acceler- 
ometers were a function of both drag and lift forces given by 


f/j pu^a - 


= I -r; si 

L 


sin Y cos X ^ 


1 - sin2 Y cos 2 x 
cos 0„ 
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•[ 

.[ 


- Cjj sin Y sin X - Cj 


sin^Y cos X sin 


cos 0, 


in X 
J 


jsin Y cos Y cos P 
- Cq cos y “ CjJ cos 0 


] 


where i, j, and k are unit vectors in the directions of the three axis of 
the accelerometers with i along the axis of the cylinder, k is in the 
direction of the nominal spin axis and j is orthogonal. Since the data 
being used in this analysis has been averaged over spin cycles, the 
component of force as expressed above were integrated over the angle X. 
The results after averaging over one spin cycle were 


Fy = siii y ^ P TTr^ (Oj, p^, y) 



cos Y 


1 _ 

2 


P U TTr C, 




Pj.Y 


where is the integrated force coefficient. These results show that 

Fy and measure the identical forces except for the factor sin y and 

cos Y* Tliis property was used by Fess and Young to obtain the angle Y 
which the spin axis malces with the velocity vector. 

Y = cos-* + 


The force coefficient Cp' v/as found by fitting a 3^^ order polynomial 


to the values of Cp and 


CjL ,2 three values 


of P; = 0, 1, and 2. 


where 


2 . 3 

Cp=A4/i-Ci:j (G+HPj + QPj +PPj ) 
A = -2 -4 Aj^ E(Y,^) /TT 
G = -4 Aj^ E (Y, J)/3 tt 


H = 4 F - 2G - 2A 


Q = -4 F + 5g/4 + 11 a/4 

P = F - g/4 - 3 a/4 

F =”4/3 - tiAj^/ 2 
n 

where E(Y, ^ complete ellcptic integral of the second kind resulting 

from the average over one spin cycle. 

The accelerometer in the x direction did not function so only Fy and 

F are treated in the analysis, 
z ^ 


18 



4. Method of analysis. 

The objective of the analysis is to find the best values ofCZj» Pj and 
density which explains the observed accelerometer output of the OVI-15. 
Although 0£j, Pj, and P are considered as unknowns in the analysis, 
certain assumptions on the characteristic variation of p are made to 
facilitate the analysis. The assumption made is that the density varia- 
tion is symmetric with respect to the perigee of the orbit. The absolute 
value of density is still treated as an unlcnown quantity. 

4.1 Least squares fit. 

Assuming symmetrical density variation about perigee, differ- 
ences in forces measured at points equal distance from perigee must 
be due to changes in the aerodynamic force coefficient, Cp. Since the 
density isr equal at these two points, we can write 

Pj (T-Atj) U2^P2(T +^^0 


where the subscript 1 indicates approach to perigee, subscript 2 indicates 
recession from perigee and T is the perigee passage time. The aero-- 
dynamic properties and forces measured at these two points must then 
satisfy the following relationship 

i 1 ^12 

Cp(0!j, Pj. Yii) Cp (eXj, Pj.YiZ^ 

where and Yi2 angles of orientation at T>»Atj^ and T +At^ 

respectively, and F =/ k/ • In the analysis, the quantity DEL^ is 

found from the proceeding relation, defined as, 


DEL- 


= i' 


' il 


i2 




where i is used to indicate a comparison made at T i Atp -A solution in 
the least squares sense is obtained by finding the values of otj and Pj 
which provide a minimum to the sum-of~DELj|^-> squared for a number of 
observations near perigee 

n 

SUM = E, (DEL-)2 
1=1 >- 


The best values of and P- are those which satisfy 

J J 


^ B(SUM) 

8(/i-sry ■ 

^ 3 

in the region of 0 Pj <. 2, and 0 ^ ^ 
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4, 2 Perigee passage time. 

The analysis requires first that the aerodynamic properties be 
different at the comparison points used in the analysis. If the Cp were 
not different, then the last equations would be satisfied for all values of 
m and Pj. Although the OVI-15 satellite was designed to maintain a Y = 

90° thoughout the orbit, considerable uncontrolled drift in the spin axis 
was found to occur. This malfunction was desirable for purposes of this 
analysis since the angle If changed considerably during a given orbit. 

This factor resulted in a changing value of Cp over the orbit which pro- 
vided a good sampling of Cp and P values for a wide range of angles of 
orientation. The analysis also requires that the perigee passage time 
be accurately known since the comparisons are made at equal points on 
each side of this time. Since the report of Fess and Young did not provide 
a perigee passage time, it was necessary to calculate that time from the 
accelerometer output. Since the aerodynamic properties are changing 
during a perigee pass due to the changing angle of orientation, the peak 
in the F curve is shifted in time from the peak in dynamic pressure. 

The derivative of F during a perigee pass is 

F ' = S Cp + s' Cp 

where S is the dynamic pressure. At the perigee passage time, the 
dynamic pressure is maximum and S = 0. Therefore, at the perigee 
passage time 

F(T) 

F (T) = Cf(T) Cj, (T) 

This equation was employed to find the true value of T for each data set 
employed. The quantities Cp and Cp are a function ofctj and Pj in 
addition to the angle Y. The analysis was able to take into consideration 
the expected shift in F output which was found to vary from 3 to 15 seconds 

depending upon CX>, p. and the rate of change in the angle ^ . 

J J 

4. 3 Speed ratio effect. 

As discussed by Karr 197 2 and Karr and Smith 197 2, changes 
in speed ratio with altitude result in a systematic increase of CD with 
altitude. The amount of increase was found to be a function of the 
satellite shape and orientation. This factor was taken into consideration 
in the analysis of data of the OVI-15 by approximating the expected change 
in aerodynamics properties with respect to speed ratio. Using information 
from Karr and Smith 197 2 and talcing into account the average over a spin 
cycle, the following speed ratio correction factor was obtained. 

COF = 1 + . 68 2/s +.. 56148/s^ 

+ . 4 sin 2 Y ( 1 . 66/s - . I528/S^) 
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The correction factor is a function of the angle of orientation which must 
be taken into account in the determination of Ofj and p.. For the altitudes 
of interest for the OVI-15, speed ratios of about 10 are obtained which 
result in an increase in the aerodynamic force coefficient of about 5%. 

The determination of density is then strongly influenced by the speed 
ratio effect. 

5. Results. 

5, 1 Data used in the analysis. 

An example of the data employed in the analysis is shown in 
Figure 4. Data from orbits number 890, 893 and 896 were employed in 
the least squares fit to Q!j and Pj. These orbits had significant changes 
in Y during perigee passage and experienced approximately the same 
atmospheric conditions. These nearly polar orbits occurred in mid- 
September 1969 with perigees at 150 km, perigee latitude at 10 ®S latitude 
with a local perigee time at about 1930. Since the sun declination 
was near +3®, the variation in density with latitude near perigee is 
approximately zero according to the Jacchia 1971 model of the atmosphere. 

The choice of orbits used in the analysis contributed to the 
reduction of errors resulting from any nonsymmetry of density variation. 
Further reduction in this type error was made by using only data within 
±10® of perigee. In addition, since data from three orbits was used, 
the errors due to wave motion or other short duration density disturbances 
would contribute only to the random error. 

The values of accelerometer output F, angle of orientation Y , 
and time in seconds from the beginning of data transmission are given 
in Table I. The units on F is in counts in which 5. 3 counts equals 10"°g 
of acceleration. The angle ^ is given in degrees. The data is seen to 
cover an angle of attack range of about 28 degrees. All the data falls 
within 150 seconds of perigee which for these orbits means that the data 
is taken within ± 10 ® of perigee. Since the true perigee is always within 
± 20 seconds of the peak F, corrected perigee times will not change the 
range of data significantly. 

5. 2 Gas Surface interaction parameters. 

Using the data given in Table I and taking into account the perigee 
passage time correction and the speed ratio correction, the results of 
the sum -of-the- squares -of- DELi are given in Figure 5. These results 
show a uni que minimum of the sum-of-the-squares-of-DELi at P j = *44 
and Y 1 -Cij = . 6. These values for 0!j and mean that the reflection is 
between a specular and diffusive type reflection in direction and has 
moderate accomadation of energy. 
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5.3 Aerodynamic properties. 

Using CXj = , 64 and P- = .44, the aerodynamic force properties 
of the OVI-15 satellite may be found using the equations already derived. 

The results are shown in Figure 6 which gives Cp as a function of angle 
of orientation. The level of Cp is dei>endent directly on the reference 
area, A, choosen to represent the satellite. The plot is given for two 
acceptable areas (1) the maximum cross sectional area seen by the flow, 
Trr^+ZrL and (2) the minimum cross sectional area seen by the flow, 

These results are for the infinite speed ratio case and must be modified 
according to speed ratio influence. 

. At angles of 0, 90 and 180°, the quantity Cp is equal to the drag 
coefficient. At all other angles, Cp is influenced by both drag and lift 
coefficients. These results show that the drag coefficient is higher than 
the value of 2.2 normally assuined in drag analysis. The speed ratio 
correction will cause these values to be inereased by about 5% to 10% for 
the altitudes at which the data was taken. 

The atmospheric density values nmy be obtained since the values 
of Cp throughout the orbit have been calculated. Assuming an orbit of 
c = . 113 and perigee at 150 km, the value of U at the data points are 
obtained and density values arc given by 

P(I, J) = 2am F{I, J)/u^(l, J)A Cp(I, J) COF 

wliere COF is the speed ratio correction factor dependent upon the angle 
Y (I> J)» 9- is the conversion factor needed to convert accelerometer counts 
into accelerometer values (a = lO^g/5.3 counts), and m is the satellite 
mass = 214 kg. The speed ratio correction requires an estimate of the 
speed ratio. For the date, time, and region of the atmosphere for which 
the data corresponds, an estimate was made of the exospheric temperature 
from information provided by Smith, 197 2. In this region of the atmosphere 
the exospheric temperature remains essentially constant and was taken to 
be 1100°. Using the Jacchia 1971 model atmosphere, a value of t/m versus 
altitude were fitted to a polynomial over the altitude range of interest from 
140 km to 220 km. The speed ratio is then given at each data point by 

S = U/J 2r t/m , 

where R is the universal gas constant. Values of density calculated in 
this manner arc about 5 to 15% less than those predicted by the Jacchia 
1971 model. A more complete discussion of these results will be made 
at a later date. 

6. Discussion of results. 

The results of this analysis are important for a number of reasons. 

First, the analysis illustrates a new method for the analysis of satellite. 
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dynamic response, Second, the results obtained forCfj and Pj are the 
most accurate of measurements of the gas surface interaction parameters 
at satellite velocity. Third, the values of Cp obtained in the analysis 
are the most accurate of measurements of satellite aerodynamic proper- 
ties, Finally, the results obtained for density are the most accurate 
of measurements of absolute values of upper atmospheric density. 

Past drag analysis have required that critical assumptions be made 
on Cp) or-p or some of the gas surface interaction parameters. The 
analysis presented here on the other hand has employed very few assump- 
tions in comparison. The assumption of symmetrical density variation 
about perigee is most subject to error. However, the possible error 
introduced by nonsymmetry is expected ‘to be much less than the errors 
committed in past drag studies. - The method presented here is far 
superior in terms of error than previous methods. 

Improvement of the errors in the present analysis could be made 
by a more accurate treatment of aerodynamics and more accurate 
measurement of accelerations and angles of orientation. The accuracy 
of angle measurements was about i 1° while the force measurements 
were accurate to i 5% for the data used. The aerodynamic description 
of the OVI-15 could be improved by employing a more accurate expression 
for the variation in -Cp' with I’j. The polynomial appr oxiination employed 
in the analysis could be improved or the exact expressions could be 
employed at the expense of computer time. 

The results obtained for Cjr as a function of angle of attack for the 
OVI-15 is of special interest because of the many drag analysis which 
have been performed on the satellite. For example, Champion, Marcos, 
and Mclsaac, 1970; Marcos and Champion, 1972; Marcos, Champion, 
and Schweinfurth, 1971; have analyzed the accelerometer data of the 
OVI-15 to reveal a number of properties of the upper atmosphere. In 
these analyses, accelerometer data was used only^ when the satellite 
was broadside into the flow. This instantaneous attitude would correspond 
exactly to the 0° or 180° spin axis orientation of the OVI-15. At this 
attitude Cjj is equal to Cp as given in Figure 5 and would have a value 
of 2.5396 for the case of infinite speed ratio and cCj = . 64 and Pj = . 44. 

This value of Cp) is 15.4% greater than the value of 2.2 which was employ- 
ed in these analysis. Additional correction would have to be made if 
speed ratio effects were. taken into account. These corrections would 
result in substantial decrease in densities reported using a Cj^ of 2. 2. 
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OVl-15 orbital data has, been analyzed by Ching 1971 and King-Hele 
and Walker 1969. The King-Hele and Walker. analysis employed a 
constant of 2.2 independent of the satellite orientation. The reference 
area used by King-Hele and Walker was midway between the maximum 
of 2. 578 nr^ and the minimum of rrr^ shown in Figure 5. On the basis 
of the reference area employed by King-Hele and Walker, a Cjq value of 
between 4.543 and 3. 889 would correspond to the values and Pj 

found in the analysis. 

The orbital decay analysis reported by Ching 1971 includes a factor 
to represent the changing aerodynamic drag properties of the OVI-15. 

The factor is based on the changing cross sectional area seen by the flow. 
The drag coefficient is considered constant while the reference area used 
in the analysis is changed by as much as 25%. A 25% correction factor 
is too large in view of the results of Figure 5 which show that the maximum 
change in Cj) A would be 16%. 

In addition, the effect these results have on past analysis of OVI-15 
data in particular, the results indicate that the assumption of Cq used 
in most drag studies have been too low. The value of Cjj = 2.0 or 2.2 
which has been used for most past drag analysis is lower than could be 
expected for most shapes with a j = .64 and Pj = .44. A sphere for 
example would have a Cjj , = 2.352 which is 7% higher than the 2 ,2 
value often used. It shouldbe noted, however, that the results obtained 
here are for one satellite surface and one atmospheric composition. It 
is expected that other surfaces and other compositions should change the 
values of aj Pj and result in a change in aerodynamic properties. 

More data must be collected before a firm value of otj and Pj can be 
assigned to a given gas and surface combination. Future work should 
be directed towards this goal. 
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TABLE I 


Orbit Number = 890 Time of Peak F = 672 


I 

F(i, 1) 

GAk4A(i, 1) 

F(i,2) 

GAMA(i, 2) 

TIME(i, 1) 

TIME(i, 2) 

1 

317.5 

137. 75 

320. 0 

134. 8 

647 

697 

2 

307. 0 

139. 50 

306.0 

133.25 

622 

722 

3 

291. 5 

141.75 

291.0 

132.25 

597 

747 

4 

270. 0 

142.50 

268.0 

130. 0 

572 

772 

5 

251.0 

144. 0 

244.0 

129. 0 

547 

797 

6 

225. 0 

145. 75 

214. 0 

127. 5 

522 

822 


I 

Orbit Number = 893 
F(i, 1) GAMA(i, 1) 

F(i,2) 

Time of Peak F = 681 
GAMA(i, 2) TIME(i, 1) 

TIME(i,2) 

1 

292. 0 

14.1. 75 

292. 0 

137.25 

650 

712 

2 

282. 0 

142. 75 

283. 5 

135. 0 

625 

7.37 

3 

264. 0 

144. 50 

263. 0 

132. 75 

600 

762 

4 

243. 0 

146. 20 

239. 5 

132.0 

575 

787 

5 

218. 0 

147. 50 

215. 0 

130. 0 

550 

812 

6 

199. 0 

148. 20 

190. 5 

128. 0 

525 

837 


I 

Orbit Number = 896 
F(i, 1) GAMA(i, 1) 

F(i,2) 

Time of Peak F = 675 
GAMA(i, 2) TIME(i, 1) 

Time(i, 2) 

1 

306. 0 

133. 5 

302. 5 

129. 5 

650 

700 

2 

302. 0 

135.0 

291. 0 

127. 7 

625 

725 

3 

289. 0 

136. 25 

275. 0 

126. 75 

600 

750 

4 

270. 5 

137.5 

256.0 

125. 0 

575 

775 

5 

250. 0 

138. 0 

231.0 

123. 0 

550 

800 

6 

226.5 

141. 0 

200.5 

120. 5 

525 

825 
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Figure 1. The gas surface interaction and the forces of the interaction. 
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CHAPTER III 


AERODYNAMIC LIFT EFFECT ON SATELLITE ORBITS 

A paper submitted to the AIAA 13th Aerospace Sciences by 
6. R. Karr, G. Cleland and L. L. DeVries. Presentation was made 
during meeting held at Pasadena, California, January 20-22, 1975. 
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AERODTHAHIC LIFT EFFECT OM SATBIXIIK OUIXS^ 

Dr. 6t ft. Kerr* •ad J. G. Clalaad** 
Dnlvaralty of ftlabna In Himtavllla 
Huntavllla. Alabau 35807 
and 

Dr. L. L. DaVrlaa*** 
NASA-Harahall Spaca Flight Cantor 
Huntavllla, Alabna 


Abatract 

Numerical quadrature la employed to obtain 
orbit perturbation reaulta from the general per* 
turbatlon equations. Both aerodynamic lift and 
drag forces are Included In the analyals of the 
satellite orbit. An exponential atmosphere with 
andjwlthout atmospheric rotation Is used, A com- 
parison Is made of the perturbations which are 
caused by atmospheric rotation with those caused 
by satellite aerodynamic effects. Results indi- 
cate that aerodynamic lift effects on the semi- 
major axis and orbit Inclination can be of the 
same order as the effects of atmosphere rotation 
depending upon the orientation of the lift vector. 
The results reveal the Importance of including 
aerodynamic lift effects In orbit perturbation 
analysis . 

I. Introduction 

The perturbations of a satellite orbit 
caused by the interaction of the satellite with 
the earth’s atmosphere has been a topic of con- 
siderable Interest since orbital flight was pro- 
posed. The forces acting on a satellite during 
Its passage through the atmosphere at speeds of 
near 8000 m/sec are predominantly the drag force 
which we will define to be that force which acts 
parallel to the velocity vector, V, of the sat- 
ellite with respect to the atmosphere. For a 
satellite having a drag coefflclant C. and a ref- 
erence area X, the aerodynamic drag force, S, is 
given by 

'd = + p ^ s ^ T?r 


where p Is the density of the atmosphere. An- 
other aerodynamic force which may act on the 
satellite Is the aerodynamic lift force which 
we define as the force perpendicular to velocity 
vector of the satellite with respect to the at- 
mosphere. Aerodynamic lift forces arise when a 
nonspherlcal satellite travels through the at- 
mosphere at an attitude such that atmospheric 
molecules are deflectad by the satellite In a 
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National Aeronautics and Space Administration, 
Marshall Spaca Flight Center, Huntsville, Ala- 
bama, through Contract NAS8-28248. 

*Assistant Research Professor, Mechanical Engi- 
neering Department, Associate Meiid>er AIAA. 
^Master's Degree Candidate, Mechanical Bngl- 
naarlng Departmant. 

***8pace Science Laboratory, NASA-MSFC, Member 
AIAA. 


nonayessetrle pattern with respect to the velocity 
vector. Consider, for exanple, a cylinder or 
cone with an axis of synmatry along the direction 
of the unit vector a. If such an object were to 
travel through the jtmosphere with a In the 
sasM direction as V, then only drag forces would 
result. However, If a were at an angle 6^ with 
respect to then a lift force will arise given 
by 

V sin 6. 


where_ C, Is the lift coefficient of the object 
and A Is assumed the same as the X used In 
the drag equation (Eq. 1). 

The dreig and lift forces given In equations 

(1) and (2) are both defined using the velocity 
with respect to the atmosphere as distinguished 
from the Inertial velocity of the satellite v. 
The Inertial velocity Is defined here as the ve- 
locity given by the orbital elements of the oscu- 
lating orbit at the satellite. 

V = (3) 

F a (1 - e cos E) ' ^ 

where a and e are the osculating elements of the 
orbit and E is the eccentric anomoly. Since the 
atmosphere at orbital altitudes may have motions 
with resgect to Inertial space, the Inertial ve- 
locity V is not necessarily the same as the ve- 
locity of the satellite with respect to the at- 
mosphere This difference Is due to the 

tlon of the atmosphere, ^ , at the satellite. 

The velocity Is then modified by the atmospher- 
ic motion to give 

= D - (4) 

Orbital perturbations under the action of 
pure drag (l.e., f a 0), with and without atmos'- 
pherlc motions with respect to Inertial reference, 
have received considerable attention for the pur- 
poses of (1) predicting the orbit of a satellite 
Into the future (dee, for example. Ref. 1) and 

(2) deducing atmospheric properties from observed 
perturbations (see, for example. Ref. X). Ibe 
Influence of nonzero lift (^ A 0) has received 
little attention due to a nuwer of factors which 
have been cited to reduce lift effects to negli- 
gible values. For example, random tunbllng of a 
nonspherlcal satellite Is often cited as a reason 
for neglecting lift since the lift vector would 
be randomly oriented and tend to average to zero 
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under such conditions. Second, lift forces ere 
thought to be negligible for satellites based on 
the argument that the aatalllte is a poor reflac- 
tor of molecules since the gas surface Interaction 
is assumed to be Inelastic. Hoirever, the assump- 
tion of Inelastic reflection Is not verified and 
the possibility still exists of significant lift 
coefficients. While many of the earlier aatal- 
lltes were nearly spherically aynoNtrlc and ran- 
domly tumbling (uncontrolled altitudes), many 
recent satellites are now spherical and may have 
large wlng-llke configurations (such as Skylab) 
and generally require altitude control by, either 
active or passive means. Such satellites suty 
then experience lift forces of greater magnitudes 
than experienced by past satellites. 

Another reason that lift* forces have re- 
ceived little attention Is the possibility that 
orbital perturbations that have been caused by 
lift forces may have been wrongly attributed . to 
other effects, such as atmospheric motions. At- 
mospheric motions that are liot collnear with the 
Inertial velocity can give rise to pure drag 
forces with components perpendicular to the in- 
ertial velocity vector. One of the major motions 
of the upper atmosphere often Included In drag 
studies. Is that movement which Is correlated 
with the rotation of the earth. The atmosphere 
Is assumed to rotate In Inertial space with 
about the same angular velocity of the earth. A 
pure drag satellite In an Inclined orbit would 
then experience components of the drag force 
which are perpendicular to the orbital plane. 

These forces cause the orbit plane to rotate In 
space. The observation of rates of change of 
the orbital plane of selected satellites has 
been used as evidence of the rate of rotation 
of the upper atmosphere under the assumption of 
pure drag only (see Ref. 3). We will show that 
aerodynamic lift forces can give rise to orbital 
plane changes of magnitudes comparable to those 
caused by rotation of the upper atmosphere. 

This Is not meant to imply that there Is no at- 
iBOspherlc rotation because there Is strong evi- 
dence to suggest a co-rotatlng atmosphere for 
the earth. We will show, however, that lift 
forces of seemingly small magnitude can at least 
produce orbital plane changes that are of the 
same magnitude as those assumed to be caused 
by atmospheric rotation. On this basis we feel 
that lift effects should be examined carefully 
In the reduction of orbit perturbation data for 
the purpose of deducing upper atmospheric mo- 
tions. 

The purpose of this paper Is to examine the 
essential characteristics of orbit perturbations 
that result when a satellite has lifting prop- 
erties. Comparisons will be made with pure drag 
cases where appropriate and comparison will also 
be made with the effect of the atmosphere ro- 
tating at the earth's rotation rate. The basic 
aerodynamic relations will be presented and an- 
ployed In the general perturbation equations. 

Lift forces will be divided Into two classes: 

(1) lift forces In the orbit plane and (2) lift 
forces perpendicular to the orbit plane. Mumar- 
Ical Integration of the perturbation equations 
was made using an exponential atmosphere for pur- 
poses of Illustration of the aerodynamic lift 
affects. Special emphasis Is given to the de- 
pendence of lift effects on high eccentricities 
In order to expand on previous work valid at low 


eccantrlclties. High eccentricity orbits also 
hava the affect of concentrating all the aaro- 
dynsnle affects at the perigee - region alnca the 
atmospheric density drops off exponentially away 
from perigee. 


II. Aerodynamic Lift 

In the free molecule envlroomant at orbital 
altitudes, aerodynamic lift and drag forces are a 
direct function of the Interaction of the atmos- 
pheric gas molecules with the exposed surfaces of 
the satellite. The characteristics of the gas 
surface Interaction for collision of upper atmos- 
pheric molecules with satellite uterlals at ve- 
locities of the order of 8000 m/aac Is not well 
understood. The gas surface Interaction Is de- 
scribed In. a^ueral manner using a modal devel- 
oped by Karr xn'^fhlch the force acting on an ele- 
ment of surface exposed to the flow Is divided 
Into two consonants. With reference to Tlgure 1, 



FIGURE 1. Illustration of Forces Acting Due 
to Gas Surface Interaction 


the force acting on the surface will have compo- 
nents F, which is the force associated with the 
momentum carried by the Incoming molecules and ' 

F* which Is the reaction force associated with 
mSlecules leaving the surface. The group velocity 
of molecules leaving the surface, D,, Is assumed 
to be some fraction of the group velocity of the 
Incoming molecules, U, due to possible Inelastic 
collisions with the surface. The relationship 
between Uj and U Is given by 

«diere_the proportionality term Is expressed as 
VI • In order that the parameter a. resemble 
what la customarily tensed the thermal nccoomo- 
datlon coefficient. A second parameter Is Intro- 
duced to provide the direction of the reflected 
force. The angle of reflection lew chosen Is a 
linear relationship between the angle of Inci- 
dence, 0, and angle of reflection, 6j, given by 

- I Pj + (1 - Pj) e (6) 

where Is the adjustable parameter with Pj - 0 
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giving specular reflection (6. =9) and F s 1 
giving diffusive reflection (0. =n/2> Independent 
of 8). The mass ‘flux Impinging on an element of 
surface Is given by 

m =-pU.n dA (7) 


where it Is the unit outward normal of the sur- 
face. The exchange of momentum that takes place 
at the surface gives the net force acting on the 
surface 


f = - (0 - Oj) p . it d A (8) 

The total drag and lift forces acting on an ob- 
ject are obtained by Integrating equation 8 over 
the exposed surface of the satellite. Detailed 
results for spheres, cylinders, cones, and flat 
plates are given In reference 4 and 5. 


uncertainties. In upper atmospheric density, accu- 
rate determinations of the gas surface Interaction 
from drag studies aloi^.ls not feasible. Recent 
work by Reiter and Mod,'\owever, on the analysis 
of drag and torque acting on a paddlewheel satel- 
lite provides the best values of gas surface In- 
teraction parameters since atmospheric density Is 
eliminated In the analysis. Reiter and Moe, em- 
ploying a number of gas surface Interaction models, 
concluded that the energy accommodation a. Is 
In the range of .75 - .95. This would putvl-a,' 

In the range of .5 to .22. The accomnodatlon 
efficients coming from their studies are high. 
Indicati ng Inelastic collisions. The value of 
Is seen to be significant even at high 
accomsMatlon c oeffici ent values (for example, 
a, = .99 gives Vi - a| = .1) which Implies that 
lift forces may also nave significant values. 

The ratio of lift forces to drag forces act- 
ing on an object Is a nondlmenslonal measure of 
the magnitude of lift effects. For the flat plate, 
this ratio Is obtained from equations 8 and 9 


The above equations show that forces perpen- 
dicular to the velocity vector can only arise 
through the TT. term In equation 8 . Since .the 
magnitude and'^dlrectlon of U are functions of 
a. and P , the lift force acting on any object is 
directly^ related to the gas surface Interaction 
parameters. This observation provides added mo- 
tivation for the study of llft-lnduced orbit per- 
turbations, since measurements of such perturba- 
tions could yield Information on the gas- surface 
Interaction. 


For purpose of Illustration of llft-lnduced 
orbit perturbations, flat plate drag and lift 
properties will be used, given by 


C 


D 


2 Bln 6 

a 


sin 0^ cosg- Pj+(2-Pj) 9^ 

(9) 

Cl = - sin 9^ Sin[-^ ^J-"<2-P^>9s] 

( 10 ) 




J sin 



P, + (2 - P.) 0 


a. 

j 


' "si 

■_n_ P + (2-P 5e. 

2 J J ® 


( 11 ) 


The lift to drag ratio given in the above equa- 
tion Is tabulated in Table I for four sets of gas 
surface Interaction parameters. The values are 
tabulated as a function of angle of attack and 
show that (1) small angles of attack of a flat 
plate yield the highest L/D values, (2) specular 
reflection causes the highest L/D values (L/D = 
tan 0 , for specular reflection) and (3) highly 
Inelaitlc collisions (a, = .99, P = 1) cause 
lift forces which are almost 10% of the drag 
force. The values of a, = .75, P = 1 were cho- 
sen to simulate the gas^surface Interaction re- 
sults obtained by Reiter and Moe. This set of 
parameters produced L/D values of near 0.5 and 
this value of L/D will be used to Illustrate the 
orbital perturbation caused by lift forces. The 
Intermediate case (a, = .5, P. = .5) was chosen 
only to Illustrate tne effect^of reflections 
other than specular and diffuse. The L/D values 
for the Intermediate case are seen to be near 
unity at the small angles of attack. 


where 0^ is the angle of attack >f the flat 
plate. For 9^ = 0 the plate Is edge-on Into the 
flow while for 9 _ tt/2, the plate has the out- 
ward surface normal In to the f low. Equation 10 
shows that the factor >/l - a.' plays an Important 
role In determining the llft^magnltude. If the 
molecules have an Inelastic collision with the 
surface, the velocity of reflection may be ex- 
tremely low, the value of a, would be near unity, 
and the lift very small. The perfectly elastic 
or perfect specular reflection (F c 0, a. = 0) 
would provide the highest possible values'^of 
lift. 


Laboratory experiments at velocities corre- 
sponding to satellite velocities are not con- 
clusive on the gas surface interaction to expect 
at orbital altitudes. Work by Hulpke^ for ex- 
ample tends to show elastic collisions while 
other experimenters hove found diffusive type 
reflection (Ref. 9) or both (Ref. 6). Due to 


A typical satellite would not have the high 
L/D values as experienced by the flat plate which 
Is an ideal lifting body at orbital altitudes. 

For comparison, Table II gives the L/D values 
obtained from results given by Sentman (Ref. 10) 
for a typical satellite shape and typical gas 
surface Interaction values. The shape Is a cy- 
linder with a conical end and having a total 
length of 4 times the diameter. The values of 
L/D are generally lower than the flat plate val- 
ues because this object maintains a high drag 
profile at all angles of attack. The peak value 
of L/D of .044 Is seen to occur at near 25*^ angle 
of attack. While this L/D Is much smaller than 
flat plate values, the lift force will be nearly 
5X of the drag force for this typical satellite 
'.shape. 


Having now established the likely existence 
of aerodynamic lift forces of magnitudes from a 
few percent of the drag force and greater, the 
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X 


( 12 ) 


ozblCal psrturbatlona that ara charaeterlatle of 
thasa forcea will ba ptaaantad. In view of tha 
f tnat tha aarodjrnanle drag forca la tl^a prln> 
elpla aourca of orbltel parturbatlona, a lift 
forca of evan a faw parcent of tha drag forca la 
axpactad to ba algnlflcant. In order to enhance 
tha lift affacta and to reduce coiaputar tlaa, an 
L/D value of .5 waa choaen for moat of tha atud- 
lea. The cooputar almulatlon employed a light 
aatelllte weight alto for tha purpoaa of enhanc- 
ing the aerodynamic affects. For these reasons, 
the results obtained ara not likely -to simulate 
any known satellite but are presented only to 
Illustrate the character of lift- Induced orbit 
perturbations. 


da 

3t 


2 f 


173 

(U a ) (1 - e^) 


d« 

It 


1 r2q-e^) cos K 
I J f r 1 - a cos B 


, 1/2 

(1-e ) sin E 


»] 


(13) 


Table I 


Flat Plate Lift to Drag Ratios 


9. 

Specular 

Diffuse 
a -.75, 
pj- 1 

Diffusa 
a, - .99 
pJ - 1 

Inter 
O,— .5 
f1-,5 



0.5 

0.1 

1.0 

2° 

28.636 

0,491 

0.100 

0,997 

5® 

11.430 

0.477 

0.099 

0.985 

K 

5.671 

0.453 

0.097 

0.947 


3.732 

0.428 

0.094 

0.896 

25° 

2.145 

0.374 

0,087 

0.772 

35° 

1.428 

0.318 

0.078 

0.642 

45° 

1.000 

0.261 

0.066 

0,514 

K 

0.700 

0.203 

0,053 

0.392 


0.466 

0.145 

0.039 

0.276 

85° 

0.088 

0,029 

a.qo6 

0,054 


Table II 

L/D for Typical Satellite Shape 
(from Sentman, Ref. 10) 


; Attack 

L/D 

0 

0.0 

10 

0.037 

15 

0.041 

25 

0.044 

35 

0.033 

45 

0.023 

55 

0,013 


III, Orbit Perturbation Equations 

The approach will be slnllar to that taken 
by others In that the expressions for perturbing 
force will ba substituted into the general per- 
turbation equations for the time derivatives of 
the osculating orbital elements a, e, 1, u), and 
fl; tha semi-major axis, eccentricity. Inclination, 
argumant of perigee, a^ right asaenslon of tha 
ascending node, respectively. Numerous good re- 
feranoas provide details of tha derivation of 
thesa aquations with specific application to drag 
forces (sea, for exmsple, Ref. 11). Tha forms of 
tha aquations most useful for this work ara given 


dl 

Tt 


' ^ :\h “ 

(n a) ( 1 - e^^) 


(14) 


iSk 

dt 


sin 

ITT 


.(... 9 . t .UI.). 


(H e)‘'‘ (1. 


•e^) 
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sin 1 


(15) 


dt 



2 sin E 
1-e cos E 


T 


(16) 



(e + cos E) N 


(17) 


r sin (9 + utf) W (18) 

1/2 1/2 
(U a) (1-e^) tan 1 


where 

r , ll/2 -1/2 

f m Kl-e^) (1 • cos E)J (1-e cos E) 

(19) 

The angle 0, the true anomaly, can be expressed In 
terms of E, the eccentric anomaly, by 


• - <”) 

and the radius, r, can also be expressed as a 
function of E by 


r - a (1 - a cos E) 


( 21 ) 


Tha T component of force Is tangent to tha orbit 
path and positive In tha direction of notion; the 
M component Is positive In the direction of tha 
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outward normal Co Che orbit path and la In the 
orbit plana; the W cooponant la the force pet^ 
pendleular to the orbit , plana with poaltlw being 
In the same direction as the orbit angular no- 
nentim vector. The aerodynamic force acting on 
the satellite will be expressed in components a- 
long », T, and W, 

Consider a flat plate In orbit with the 
orientation of the plate held constant with re- 
spect to the N, T, V coordinate system. The unit 
vector X will be taken as the normal to the flat, 
plate with components In Che N, T, W directions 
given by 


cos 0^ cos 


Sjj ■= cos sin 0^ 


NT 


( 22 ) 


e elm •hi ^(1 + e eoa 6) 

cos e aln oj 


+ cos 


a^ » |cos 0JJ eoa *HT (1 + e cos 6) 

- cos sin 0JJJ e sin oj 


( 26 ) 




= sin a 


Equations 1 and 2 give the lift and drag force as 
a function of the relative velocity V. The re- 
lative velocity ^ Is first found with respect to 
the R, S, H coordinate system In which the R di- 
rection is radially outward from Che geocenter 
and S Is perpendicular to R and W with positive 
being In the direction of motion. An atsiosphere 
rotating at the earth's rotation rate, ^ , Is to 
be considered the only atmospheric motloS for 
purposes of this study. The velocity flald pro- 
duced by the atmospheric rotation Is given by 


The angle 6^ called for in the drag and lift re- 
lationships can be obtained In terms of the orbit 

parameter and the anglea 0 and 0 from 

W NT 

®» ^ i vr*) <27) 


where 


IT - r R + (r 0-n^ r cos 1)S +[n^ r cos(e 4 w) 

sin 1 w] (28) 




= r n [cos IS- sin 1 cos(04w)w] (23) 


It Is convenient to write the velocity ^ as a 
function of E and E. To do this, the following 
relationships are employed 


where S, W, and r are unit vectors. The Inertial 
velocity of the satellite Is given by 


” A -A 

T^*rR + r0S 


(24) 


The relative velocity V Is by equation 4. 

Equations 4, 23, and 24 are substituted Into 
equations 1 and 2 to obtain the components of lift 
and drag In the R, S, V coordinate system. The 
unit vector will have components a_, me, and 
^ In the R, S, W system. The transforknaMon 
from the N, T, W system to the R, S, V system Is 
obtained from Che following relationship between 
the unit vectors 


A 

N 


A 

T 


1 


[(1 + e cos 0) R - e aln 0 s] 


(25) 


r ^ 

|_e sin 6 R ^ (1 4 - e cos 


9) S] 


r 

> a e aln E E 



2 1/2 


r e 

- a (1 - e^) 

E 


5/2 


r 

" (1 - • 

cos 


( 29 ) 


.2 1 






a^^2^1 _ g go, j5j 


Substitution of these relationships Into the force 
equations, and taking components in the i, i, T 
directions, we obtain the following com p o n ents of 
force &, S, and T. 


R - p 


V^E 


E I 

[ n sin 0^ 


( 1 - e eoa I) 


Using aquations 25 and 22, we find 


sin R 


(BLCOte, +Bp)J 


( 30 ) 
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S -p 


<1- 

sin 6^ 


[ , 1/2 n cos 1 
(l.«Z) . -a-j (i,« CO# B )^J 

(-r) ®. + v} 

« - P ^ ® { n!iA, (1 - • ^o* « 

^ ■ <• 


(31) 


( 1-0 COB E) cos (8 + (u)sln 1 

(32) 


n V 

<®L ®s + V I 


C D A C A 

®D ■ ~ 2" ~ • h *= "TT~ 


of sstelllte and n = / s'* . 

This concludes Che development of the orbit 
perturbation relationships. Equations 30, 31, 
and 32 along with the transformation relations 
given In equation 25, provide the necessary re- 
lations to write the orbit perturbation equations 
In Che form 


dE “ *a^®L* ®D* '^W* **HT’ “* *’ 

S ®D* ^W* ®NT’ “• *’ 

etc. (33) 


The step size on AE was taken as small as 0.1 de- 
gree depending on the rate of variation of the or- 
bital elements. The accuracy of the results Is 
not expected to be good near the final stages of 
decay since the orbital elasients are changing rap- 
idly. Errors also build up at low eccentricities 
due to the high rate of change In (U which Invali- 
dates the assumption of constant orbital angular 
momentuB Implied In equations 29. 

The Integration of the equations stepwise 
in E Is a departure from previous work by Qook 
(Ref. 14) In which a series exptmslon was employ- 
ed to facilitate the Integration analytically 
from 0 to 2rr. We did not Cake this approach for 
two reasons; First, the results obtained by Cook 
axe only valid at low eccentricities (less than 
.2) and expansions for high eccentricities were 
not found. Second, we were Interested In Che lift 
effects on orbit decay which are not apparent If 
one Cakes, as Cook did, equation 12 to Imply that 
lift forces have no effect on decay of the seml- 
s;aJor axis. Equation 12 shows that the rate of 
change In a Is, Co first order, a function of 
only, forces tangent Co the orbital path. For a 
non-rotating atmosphere, drag forces only would 
produce tangential components and lift forces 
would not effect a. We find that In a non-ro- 
tatlng atmosphere, n =0, lift forces do affect 
the value of a ovef Chat of drag alone. This 
Is clearly a second order effect which we find Is 
comparable to the effect of atmospheric rotation 
on a for L/D = .5. Results such as these are 
found’ by stepwise Integration of the orbit equa- 
tions and would be lost If Integration from E = 0 
Co 2rr were done assuming all parameters remain 
constant over the Interval. 

V. Results 


The results are grouped Into three classes 
of lift orientation: (A) Lift forces perpendicu- 

lar Co the orbit plane, (B) Lift forces in the 
orbit plane and In Che same direction around the 
orbit, and (C) Lift forces In the orbit plane 
which change dlirectlon at perigee and apogee. 

The first and last type of lift orientation serves 
to approximate cases In which satellites maintain 
an orientation with respect to Inertial space. 

The second type of lift orientation serves to ap- 
proximate Che earth oriented object. 


The atmosphere density called for In the pertur- 
bation equations Is given by 

p (r) = pp e - (34) 

where p Is the density at perigee, h Is Che 
height of perigee, H Is the scale heigRt at 
perigee, and h Is the height of the satellite 
given by h > r - R- where R^ Is the radius of 
the earth. “ “ 


IV. Hethod of Calculations 

Hie orbit perturbation equations were inte- 
grated numerically using E as the Independent 
variable. A modified Runge-Kutta fourth order 
process 'due to Glll^ana discussed In ref.l3was used. 
Double precision was used to maintain accuracy. 


A. Lift Perpendicular to the Orbit Plane 

For a flat plate with angles of orientation 
with respect Co the N, T, W coordinate system 
given by (0„ = constant other than zero, = 0), 
lift forces will be produced In the W direction 
which Is perpendicular to the orbit plane. The 
effect on the orbit of aerodynamic forces perpen- 
dicular Co Che orbit plane have been studied by 
Cook and Pllomer (Ref. 15) and Cook (Ref, 16) for 
the case of a satellite having pure drag with the 
perpendicular forces arising due to the rotation 
of Che atmosphere. For comparison, we took the 
case of a non- rotating atmosphere with an aero- 
dynamic lift force of L/D .5 and .1 acting per- 
pendicular to the orbit. 


The orbit parameters most sensitive to forces 
perpendicular to the orbit plane Is the orbital 
Inclination 1 and right ascension of the ascend- 
ing node Cl As pointed out In refemnee 15, the 
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pextuxbctions In Q cnuaed by Mrodynatic foxcu 
perpendteulnx to thn ozblt «r« th«n .IX of 
those eeused by the obletenesa of the eerllh. The 
Inclination 1 la not affected |>y earth oblate- 
ness to such an extent and Is thanfore a sensi- 
tive indicator of aarodynsnlc forces perpendicu- 
lar to the plane. 

The characteristics of orbital inclination 
changes caused by aerodynanic lift were found to 
be not nueh different fcoo those caused by atmo- 
spheric rotation acting on a pure drag satellite. 
One difference la, of course, that the ma^itiide 
of component of drag foxes in the W direction 
caused by atmospheric rotation is proportional to 
ain i whereas aerodynamic lift forces are inde- 
pendent of 1. Therefore, inclination changes of 
near zero inclination orbits must be caused by 
aerodynamic lift effects only, since atmospheric 
rotation can have no Influence at zero inclina- 
tion. For comparison purposes, Che lift Induced 
rate of change in inclination with no atmospheric 
rotation was normalized with respect Co the in- 
clination rateg of change for a pure drag aatel- 
llCa at 1 B 45 in an atmosphere rotating at the 
earth's rotation rate, niese results are pre- 
sented in Flgura 2 for L/D = .5 md L/D - .1 
with Bj, » .2 m^/.Kg and = .1 mVKg. 



FIGURE 2. OuC-of-Flane Effects On Orbital 
Inclination 


The notation is used to signify that the lift 
force is in the W direction. Figure 2 shows that 
the rate of change in inclination for L/D = .1 is 
over twice that that would be caused by atmo- 
spheric rotation. The curve for L/D = ,5 shows 
Chat Che dffect is proportional to L/D. Figure 
2 also shows a comparison of the atmospheric ro- 
tation at i • 90° with respect to that at 1 ° 45°. 
The effect is proportional to the ratio of the 
sines of the angles as would bo expected. 

Figure 2 shows an increase in the rate ratio 
as a function of e. The reason for this increase 
is due Co the decrease in di/dt for the atsw- 
spherlc rotation effect as e is increased since 
Che inertial velocity near perigee increases with 


e. The atmospheric rotation effect on i de- 
creases since it is proportional to the ratio of 
V to V which decreases with eccentricity. The 
atrodynsnlc lift force is not affected by e since 
L/D is independent of velocity. 


B. In-Plane-Lift Forces Constant in Direction 


Consider new a flat place with angles of 
orientation such that fly « 0 snd s constant 
(other than zero) , This orientation would pro- 
duce a lift force that is in the N direction and 
constant in sign around the orbit. Two cases will 
be considered, (1) L/D >= .5 with the lift force 
directed always in the positive N direction and 
(2) L/D = ,5 with the lift force directed always 
in the negative M direction. As pointed out above, 
aquation 12 shows that these lift forces do not 
cause a change in a to first order. 


We investigated the change in a and e (A a/ 

A C and A e/A c) at a point near perigee, but al- 
ways after perigee, using the numerical methods 
described. The results are given in figures 3, 4, 
and 5 in which the effect of lift and the effect 
of atmospheric rotation are both compared to the 
effects caused by pure drag only with no atmospher- 
ic rotation. The notation L, is used to signify 
that the lift force is in the N direction. 


Figure 3 shows that A a/A t la positive or 
negative depending upon the orientation of the 
lift vector. The magnitude of the effect is de- 
pendent on the eccentricity, being larger for 
smaller values of e. Also plotted in Figure 3 
is the effect gn A a/A t due to a rotating atmo- 
sphere (1 e 45°) , The magnitudes of the two ef- 
fects are comparable and are seen to be of order 
5X of the drag-nonrotatlng-atmosphere effect a- 
lone. The lift effects are second order effects 
as pointed out above while the atmospheric ro- 
tation effects are first and second order. The 
first order effect due to atmospheric rotation 
comes in because the force in the T direction ^s 
decreased by the wind for 1 between 0 and 90°. 
Figure 4 shows that the effect on e has much 
the same character as the effects on a. 



InitW CeetMrtcity, • 


FZODRB 3, In-Plana-Lift Effects on Seml-Hajor 
Axis Decay Rates 
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FUnnUt 4. In-Flane-Llft Effects on Eccentric* 
Ity Decey Rates 


Figure 5 shove the dependence of 4 e/4 t on 
perigee height for two eccentricity values e > .2 
and e ■ .8. Below a perigee altitude of 130 Kib, 
the rotating atanosphere effect Is seen to tend to 
reverse while the lift effects become larger and 
do not change sign as does the atmospheric rotation 



FIGURE 5. In-Flane-Llft Effects on Instanta- 
neous Semi- Major Axis Decay Rates 
as a Function of Perigee Height 
The effect of constant lift on the orbital 
llfetlaw is shown In Figure 6. The lifetime Is 
ccoputed from time of Initial conditions until the 
altitude becomes circular. Although the computation 
technique was not deslmed to handle zero eccentric- 
ities, little error Is Involved since Che orbit 
eccentricity usually would become zero from Ini- 
tially high values only very near final decay. The 
llfetlsMS are normalized with respect to the life- 
time of a pure drag satellite having the sane drag 
as the lifting satellite. The atmosphere was taken 
to have zero rotation In order to clearly show the 
effect of lift. The results show that In-plane lift 
effects tend to decrease the lifetime for high ec- 
centricity orbits, independent of the direction of 



FIGURE 6. In-Plane-Llft Effects on Satellite 
Lifetimes 


Che lift vector N (assumed to be constant In sign 
around the orbit). AC low eccentricities, the ef- 
fect on lifetime Is less but there Is separation of 
the effect according Co the sign and magnitude of 
Che lift vector. The results at low eccentricities 
contain some error due to the rapid rotation of the 
orbit In space. For this reason, the magnitudes of 
the influence on lifetime may be In error at low 
eccenCriclCles by an unknown amount. The curve Is 
presented here In order to show the clear tendency 
of lift effects to reduce lifetime as eccentricity 
Is Increased. 

With reference to equation 16, Che effect of 
constant lift force around the orbit should have a 
significant effect on J) since the M component of 
force Is multiplied by cos E. This was also noted 
and evaluated by Cook (Ref. l4) for low eccentrici- 
ties. Using the same aerodynamic values used by 
Cook, we evaluated iL for low and high eccentrici- 
ties and compared the results as seen in Figure 7, 

The comparison with the Cook result Is seen to be 
good at the lowest eccentricities used In our numert 
leal work. In addition, our results show that <I> 
for higher eccentricities continues to decrease with 
Increasing e. 

C. In-Plane-Llft Forces Which Change In Sign at 
PerlKee and Apogee 

Consider a flat plate flown such that the 
angle of attack in the NT plane changes dlseontlnu- 
ously In sign, but remains constant in magnitude, at 
perigee and apogee. This special case is of interest 
since it has the effect of appro«lmatlng a space- 
craft with attitude controlled with respect to In- 
ertial space. Two eases arise; (1) lift in positive 
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Ecctniricily, • 

FIGURE 7. Argument of Perigee Rate of Change 
with Lift and Drag 


direction of N before perigee and in negative di- 
rection after perigee and (2) lift In negative di- 
rection of N before perigee and in positive di- 
rection after perigee. These types of lift his- 
tories have the effect of making the N sin E term 
in the equation for de/dt (Eq. 13)elther always 
negative, case 1, or always positive, case 2. 

Case 1 Is the same as that considered by Cook in 
reference 14. Instantaneous values of A a/A t 
and A a/A t at a point Just past perigee for both 
positive and negative values of lift are given In 
figures 3 and 4. Values of A e/T where T is 
the orbit period were computed and found to agree 
well with the results published by Cook and will 
not be presented here. Of special Intetest, how- 
ever, was the effect of the discontinuous lift 
cases on orbit lifetime which was not considered 
by Cook. 

The lifetime of the orbit was found to be 
strongly affected by the discontinuous lift cases. 
The results are given In Table III In which the 
lifetime Is nonuallzed to the lifetime of a satel- 
lite having drag of the same magnitude as the lift- 
ing satellite and with no atmospheric rotation. 

The lifetimes for' case 1 lift histories were seen 
to be Increased by 1.7 to 3.5 times over the drag 
only lifetimes. Only three lifetimes were com- 
puted due to the long computing times required as 

Table III 

Discontinuous Lift Effect on Lifetime (L/D=.5) 


CLlfetlme) / (Pure Drag Lifetime) 


e 

Case 1 

Case 2 


N sin E Is Neg. 

N sin E is Fos. 

.1 

1.7 

.52 

.2 

3.2 

.48 

.3 

3,5 

.41 


the eccentricity was Increased. The case 2 lift 
history is seen to reduce the lifetime to about 
half compared to the drag-only case. Careful anal- 
ysis of our results shows that the case 1 lift 


history has the effect of reducing the perigee de- . 
cay rate of the orbit, thereby reducing the drag 
effect, whereas the case 2 lift history causes the 
perigee height to decrease faster than the pure 
drag case. This effect is best understood by con- 
sidering the relationship for perigee radius, r , 
given by 

Tp = (1 - e) a (35) 


The time derivative of r is 

P 


dr 

^ 


da 

dt 


de 

dt 


(36) 


The values of da/dt and de/dt are negative for pure 
drag and for both case 1 and case 2 lift histories . 
The action of case 1 and case 2 lift histories is 
to modify the decay of perigee height with respect 
to the pure drag case in the following manner; 


dr / \ 

= (1 + X) ( 1 + (higher order terms) 

(37) 

where X Is a number of magnitude .1. The plus 
sign Is taken for case 1 and the negative sign Is 
taken for case 2. The subscript o refers to the 
pure drag values of . 

VI. Conclusions 


Satellite lift forces were shown to be of the 
order of a few percent of the drag force for an 
ordinary satellite and may range up to near unity 
for a satellite designed to have high lifting forces. 
Reference to figure 2 concerning the lift effects 
on the orbit Inclination leads us to conclude that 
an L/D of only .01 would produce an orbital in- 
clination rate of change of a magnitude nearly 25% 
of that attributed to atmospheric rotation. Clear- 
ly, then, an attitude stabilized satellite with a 
lift force consistently perpendicular to the orbit 
plane could produce large changes In the Inclina- 
tion which may be wrongly attributed to high ve- 
locity upper atmospheric winds . 

The results obtained concerning the effect 
of different lift histories on the lifetime of the 
orbit are conclusive on three points: (1) constant 

lift directed either positive or negative along N 
about the orbit causes a reduction in lifetime at 
high eccentricities, (2) a discontinuous lift with 
change In sign at perigee and apogee, such as to 
have positive lift before perigee, causes Increased 
lifetime over a drag-only satellite, (3) a discon- 
tinuous lift change of the other type. In which 
the lift Is negative before perigee, causes a re- 
duced lifetime of the satellite. These conclusions 
are apparent from the results given In figure 6 
and those given In Table III. 

Due to difficulties we had In the numerical 
procedure near final decay of the satellite, the 
numbers on lifetime may be In error. For this 
reason, we stress that the results merely provide 
the character of the orbital perturbations due to 
lift and may be In error In absolute value . We 
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found It difficult to put • reallitlc oatliMto of 
•rror on tha raaulti obtained. All wocfc waa eaiv 
tied out In double praclalon and tha atap alca waa 
dacraaaad until no changaa In pravloua raaulta ware 
found. We did find variations fron tha amodth 
cutvea of flgurea 3, 4, and 5 of tha order of 1 to 
S%. The raaulta are conaldarad accurate enough to 
Indicate tha control 1 Ability of satellite orbits by 
lift. 

The results show the influence of satellite 
lift on the orbital elements and also open tha 
possibility of utilizing satellite lift to, non- 
propulslvely control the satellite orbit plane and 
orbital lifetime. A flying spacecraft would have 
an advantage over the propulsive satellite If the 
aerodynamic shape could be Incorporated Into the 
design without Increased weight. The weight of the 
propulsion system and Its fuel- could then be used 
as payload. 

Future work should be done on the observation 
of lifting effects on existing satellite orbits In 
order to confirm the results obtained here. In 
addition, analytic work and further numerical 
studies need to be done on the Influence of lift 
on the orbital elements. 
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DUAL FALLING SPHERE DETERMINATION OF DENSITY AND 
TRANSITION FLOW PARAMETER 

Paper submitted to AIAA 12th Aerospace Sciences Meeting by 
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D. C., January 30-February 1, 1974. 
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Abstract 

A new approach to the analysis of falling 
sphere drag data Is described In which Che data 
from two trajectories through the same region of 
the atmosphere are. analyzed simultaneously. The 
analysis provides Important aerodynamic Infor- 
mation vdilch la used to obtain an Improved value 
of atmospheric density. The technique Is applied 
Co a sec of falling sphere', data In which a sphere 
transition-flow parameter and atmospheric density 
results are obtained In the 80-120 km region from 
published data for falling- spheres over Kwajaleln. 
Another set of data for a falling sphere test over 
Wallops Island is also analyzed with comparable 
results. 

Introduction 

The falling sphere technique has been the 
prime source of measurements of atmospheric 
density and temperature In the Important altitude 
range of 80 to 120 km. In particular, the results 
from three falling sphere experimental groups 
provided the Information used to supplement Che 
1962 U. S. Standard Atmosphere.**^ These groups 
were Peterson, - Hansen, McWatters and Bonfantl,’ of 
the University of Michigan; Falre and Champion'* of 
the Air Force Cambridge Research Laboratories; and 
Pearson’ of Australian Weapons Research Establish- 
ment. The results obtained by these groups are 
summarized In Che 1966 supplements to Che U. S. 
Standard Atmosphere.^ The method of analysis In 
all these experiments was to first measure the 
acceleration, a, acting on the sphere from either 
trajectory analysis or accelerometer readings. The 
drag equation 

Drag - y p V^Cd A - ma (1) 

Is than employed to deduce the density, p, where V 
Is the velocity, Cq la the drag coefficient, A Is 
a reference area and m Is the mass. The mass, 
acceleration, velocity, and area are all measured 
quantities, leaving only p and Cd as unknowns In 
equation (1) . A Cable of Cd as a function of 
Reynolds number, R^, and Mach number, M, Is usually 
employed to obtain a value of Cd- A density value 
Is chan determined by solving equation (1). Since 
Cd Is a function of Reynolds number and Mach 
number which are density and temperature dependent. 


*Thls work was supported by the National Aero- 
nautics and Space Administration under Contract 
NAS8-28248 through NASA/MSFC, Huntsville, 

Alabasw. 


-the analysis of the drag data usually Involves a 
number of Interatlons until there Is convergence 
to a final density value. 

The work reported here Is the result of an 
examination of the falling sphere technique for 
the purpose of proposing possible Improvements, 
particularly In the area of aerodynamics. As with 
all .drag deduced density experiments, , the values 
employed for the drag coefficient In the data 
reduction Is a major source of error. The possi- 
bility of error due to drag coefficient Is largest 
In the 80 to 120 km region, due to the passage of 
Che sphere through various aerodynamic regimes. A 
typical falling sphere trajectory will be In a 
free molecule flow regime at high altitude and will 
pass through transition flow Into conclnum flow at 
low altitudes during descent. In addition, the 
speed of the falling sphere may pass from subsonic 
to supersonic and back to subsonic during a typical 
trajectory. 

A considerable Improvement In knowledge of 
sphere drag coefficients has been made recently 
through a comprehensive experimental program at 
ARO which was sponsored by AFCRL for specific 
application to the falling sphere program. The ARO 
work reported by Bailey and Hiatt’ covers a 
velocity range from 0.1 to 6.0 in Mach number and 
a Reynolds number range from 20 to 100,000. While 
this data is extremely useful, there Is still a 
lack of accurate Information for the near-free 
molecule drag coefficients which would correspond 
to Reynolds numbers below 100. 

In the wcrrk reported here, two new methods 
of falling sphere data analysis are proposed 
which should find application In the low density, 
high altitude region of the atmosphere where 
accurate aerodynamic Information Is still lacking. 
The proposed methods Involve the simultaneous 
analysis of the trajectory of two spheres 
travelling through the same region of the atmo- 
sphere with different velocity. As an Illus- 
tration, the falling sphere data of Peterson,’ 
et al. and Peterson and McWatters’ Is analyzed 
using one of the proposed techniques. 

Basic Theory of Dual Falling 
Sphere Experiments 

Consider an experiment In which two sphere 
drag measurements are performed In the same 
region of the atmosphere but at different 
velocities. 
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i p V *1 <2) 

°2 “ 7 P V ^2 

Since Che spheres are In the sane region of the 
atmosphere, the values of density are the same for 
equation (2) and (3), Assume also that the 
dependence of Cq on p, V and a third quantity B, 
which will be discussed later. Is well understood. 
Then, In general 

Cdi (P. B. Vi) Cdj (P. B, Vj) (4) 

and equations (2) and (3) represent a set of two 
equations In the unknowns p and B. The values of 
p and B are found by solving equations (2) and (3) 
simultaneously. 


p-p (Di, 

^2# ^2* ^1* ^2^ 

(5) 

B - B (D^, 

°2' *1' ^2^ 

(6) 


In the process, not only Is a value of density 
obtained but also the two drag coefficients are 
determined for the two cases. 

The basic theory of dual falling sphere 
analysis depends on the fact that Che drag 
coefficient Is not Independent of the velocity and 
can be written as the function of two other 
parameters at most. The dependence of the drag 
coefficient on density and the quantity B must be 
known In order to write equations (S) and (6) . The 
quantity B will be seen to be the temperature 
for the first case discussed below and the 
transition flow parameter for the second case. 

If the drag coefficient were a function of more 
than. two unknowns, for example 

Cd - Cjj (p, B^, B^, B^, • • • Ba) (7) 

then two possibilities exist. One could consider 
multiple falling sphere experiments designed so 
that m + 1 drag measurement would provide Infor- 
mation necessary Co Invert the m + 1 drag 
equations, giving 

p-p (Di,D 2 ••• Da^.l,yl, V2 ,*”-V^jAi,A2*-*Aj^i) 
®1 " (Dj, D2, ••• Vj ) 


leas than In single falling sphere analysis where 
all Che quantities must be assum^. The approach 
taken In the following discussion Is to reduce the 
unknowns In Che drag coefficient to the density and 
one other quantity called B In the above. The 
potential of the dual' falling sphere technique 
appears to be greatest In the high altitude region 
In which Che flow Is free molecule and the drag 
coefficient is very sensitive to the atmospheric 
temperature as discussed In the next section. 

Temperature and Density 
Determination In Free Molecule Flow 


In the free molecule flow regime the drag 
coefficients of a sphere Is strongly dependent 
upon the speed ratio for speed ratios of order 
unity and less. The speed ratio Is defined as the 
ratio of the sphere velocity to the thermal 
velocity of the gas given by 

S - V/V2 RT (9) 

the free molecule sphere drag coefficient Is 
obtained from free molecule theory to be 



where K is a factor of order unity dependent on 
the gas surface Interaction. 

The drag coefficient Is found to be Inde- 
pendent of density in a free molecule flow. The 
nature of the dependence of on S Is better 
Illustrated by the expansion of equation (10) for 
the cases of large and small values of S. The 
results are 


CDf + S--^ S’l (1+K) 

TT L3 S 15 240 J 

(11) 

S < 1.25 

( 12 ) 

S > 1.25 


B2 " B2 (D]^, D2, ••• I>iB+l« V^ ) 

a 

Bm “ B„ (Di, •••• D^i, V^, ) 


(8) Equation (11) shows chat tends to Infinity 

as S tends to zero while equation (12) shows Chat 
CDfn fends to 2 (l + K) as S tends to Infinity. 
Values of S of order unity and less are seen to 
cause the greatest variation In Cg£^. 


The second possibility, and the one more 
practical Co consider, would be to assume all but 
two of the quantities needed to determine the drag 
coefficient In equation (7). An error Is of course 
made depending upon the accuracy of Che 
assumptions but the total error Is expected to be 


Since falling sphere velocities at high 
altitudes typically correspond to speed ratios 
of unity or less, a dual falling sphere analysis 
may be feasible and fruitfully applied In this 
region. Consider two sphere drag measurements 
at free molecule conditions 
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-1 -l/*! - I P V (h- h> 

- P Vi* Cdj^ (T. V^, K^) (13) 

■2 *2^2 “ -j P^2* Cdj^ (S 2, *2^ 

■ f pV <T. V 2 , K 2 ) (14) 


where the speed retlo has been written as e 
function of tesiperecure which mist be the saae for 
both aeasurenents. Since free aolecule dreg 
coefficient Is Independent of density, the density 
can be ellalneted In the above equations giving 


"1 


1 

2 


‘"fn 


ll 

(T. Vi, Kjl) Aj 


(15) 


|v*Cn (T, Vj, K 2 ) Aj 

Equation (15) now contains the unknowns T, Kj^, and 
K 2 > If the two spheres have the same surface 
properties and V^ls not much different from V 2 , one 
would expect the gas surface Interaction to be the 
same for both spheres. Therefore, 


however. Is lndepen()ent of the density determi- 
nation and becomes more accurate at higher 
altitudes where free molecule conditions prevail. 

A survey of published falling sphere data did 
not produce data of the type needed for an example 
analysis of the free molecule type. The proper 
data could be obtained by launching two spheres at 
near the same time but with different velocities. 
Another technique would be to track the sphere 
both during ascent and descent since the veloci- 
ties would be different due to drag effects. Data 
of Che latter type Is available but only at lower 
altitudes where the flow is transition rather 
than free molecule. Analysis In the transition 
regime Is considered In the following section. 

Dual Falling Sphere Analysis 
In the Transition Regime 

As a falling sphere passes Into regions of 
greater density, the drag coefficient must change 
from a free molecule value (2 and greater) to a 
continuum value (1 and less). The flow regime 
between the limits of free molecule and continuum 
Is termed the transition flow regime. Since no 
theoretical expression Is available which can be 
accurately applied to the transition flow regime, 
empirical and seml-emplrlcal relationships are 
commonly employed. One such relationship given by 
Matting’ has application to drag coefficient 
determination In the near-free molecule side of the 
transition regime. The expression Is given as 


Ki - K2 - K 


(16) Cd - CDc + (CDfn - Cpc) e"E/Kn (I 9 ) 


and Che common factor (1-t-K) can be eliminated 
from equation (IS) , leaving a single equation In 
the unknown temperature T. Therefore, from 
equation (15) 

T * T (mj^ aj^, m 2 S 2 , Vj^, V 2 , Aj^, A 2 ) (17) 

and either equation (13) or (14) may be used to 
obtain 

p - p j^m a, (T, V, K)j (18) 

Notice that the value of K becomes Important In 
determining the value of density but la not 
required In determining temperature. 


where Cp|, Is the continuum drag coefficient, Kn Is 
the Kundsen number, and E Is the parameter which 
must be determined from experiment. The above 
equation Is seml-emplrlcal based on a first 
collision analysis of near-free molecule flow. 

For application to falling sphere analysis 
the quantity E/Kn can be expressed as a function 
of density since 


where r Is the sphere radius. Therefore, write 

E/Kn H C 3 p r (20) 


Discussion of Proposed Free 
Molecule Analysis' 

The above procedure would provide a more 
accurate measurement of temperature at high 
altitudes than chat provided by single sphere 
experiments. In single falling sphere experiments, 
the temperature Is deduced from the density values 
by Integration of the hydrostatic equation 
beginning at the high altitudes. As discussed by 
Barcman, Chang, Jones and Llu,^ this method of 
determining temperature Is subject to large errors 
at ths high altitudes. The dual sphere analysis. 


where C 3 will be termed the transition flow onset 
parameter. Since and Cd,. are functions of 

velocity, temperature, and the gas surface inter- 
action, the functional dependence of Che 
transition flow drag coefficient is written from 
equations (19) and (2P) 

% - Cd (P, C 3 , T, V, K) (21) 

Since C 0 contains four unknown parameters, 
two parameters must be assumed In order to perform 
dual ttiUef sphere analysis In the transition (falling) 
..egime. The value of temperature will be assumed 
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to be given by the usual methods of Integration of 
the hydrostatic equation since Bartman,' et al. 
point out that the temperature deduced by the 
Integration methods becomes more accurate at the 
lover altitudes. The second parameter to be 
assumed In the analysis Is the gas surface Interr- 
action parameter, K. The value of K « 0 will be 
employed In the analysis for reasons discussed 
later. A gas surface Interaction In which 
molecules are reflected In the specular direction 
corresponds to K - 0, Independent of the degree 
of acconmodatlon to the surface temperature. 
Therefore, C3, Is chosen as the unknown parameter 
which will be determined In the analysis In 
addition to the density p. The value of C3 la not 
well established due to a lack of experimental 
results In the near-free molecule regime. The 
analysis will then help establish the value of 
this Important parameter for use In future experi- 
ments. 

Having chosen the unknown parameters, the 
method of analysis Is similar to the free molecule 
flow analysis. Consider two drag measurements 
In the transition regime at the same region of the 
atmosphere but at different velocities. 

Di - I p Vi* |'0 d^j + 

“2 " i P '^2 [^Dc 2 ''' 

where the spheres are taken to be the same size so 
that Aj “ A2 “ A and r^^ “ t2 “ r. The quantity 
exp (-C3Pr) can be eliminated In the above 
equation giving a single equation In the unknown 
density p. 

(Di/^V*A) (CD f„2-CDc2)-<“2/|v|A) (CDf^i-Cp^l) 

p — = (23) 

^°cl ^Dfm2" ‘^°c2 '^Dfml 

and C3 Is found from either of equations (22) once 
a value of p Is determined from equation (23). 


)e-C3PrjA 

,( 22 ) 

<%m2-%2 )«"^3P>jA 


C3 



CD /I A) ^ - Cp^ 


(24) 


The transition flow analysis has been 
successfully applied to five seta of falling sphere 
measurements over Kwajaleln made by the University 
of Michigan group In 1963 and 1964 as described in 
the next section. 


Analysis of Kwalaleln Falling 
Sphere Measurements 

The University of Michigan falling sphere 
measurements consist of data taken both during 
ascent and during descent for one of the three 
spheres ejected from the rocket during the ascent 
phase. The spheres were made of Mylar Inflated 
with Isopentane having an Inflated diameter of 


0.66m and a mass of 50 grams. The velocity 
altitude history of a typical flight Is shown In 
Figure 1. In this particular flight there la an 
overlap of ascent and descent data In the region 
between 90 and 110 km. Data of the type shown In 
Figure 1 are suitable for analysis as a dual fall- 
ing sphere experiment using the method outlined 
above . 


Unfortunately, not all flights had regions of 
overlap. Of the 13 successful flights over 
Kwajaleln, only six have any overlap. The sounding 
number and the region of overlap of the ascent and 
descent are the following: 


Sounding 2; 
Sounding 3; 
Sounding 8; 
Sounding 12; 
Sounding 13; 
Sounding 14; 


100 to 102 km 
99 to 102 km 
104 to 109 km 
99 to 104 km 
96 to 107 km 
at 102 km 


The plots of these data are given In the reference 
by Peterson,’ et al. The data used for the work 
reported here was obtained In the form of computer 
output of the University of Michigan analysis, a 
sample of which Is also given in reference 3. The 
temperature data obtained In that analysis was 
used directly In this work after smoothing over the 
data In the region from 80 to 120 km using a third 
order least squares polynomial. The drag term 
required In equations (23) and (24) was obtained 
using the published values of Cjj and p in the 
following relationship 


D 



A 


(P Cp) 


published values 


(25) 


Due to noise In the data, the application of 
equations (23) and (24) could not be made point 
by point. A smoothing of the data was then 
performed over the region of 80 to 120 km employ- 
ing a second order least squares polynomial fit. 
This meant that the ascent and descent 
trajectories were smoothed over a distance of about 
20 km. 

The Cq, values used were obtained from either 
equations (11) or (12) and Cd^, values were obtained 
using a polynomial fit to the data given In 
Reference 6 for the highest Reynolds numbers. The 
mean molecular weight needed to determine speed 
ratio or Mach number was obtained from a polynomial 
fit to the mean molecular weight given In the 1962 
U. S. Standard’ for the region between 80 and 120 
km. 


Discussion of Results for 66 am Spheres 

Values of were calculated In the region 
of overlap for each of the soundings. Due to the 
smoothing operations only one value of Cj could 
be obtained from each sounding. Soundings number 
8, 13, and 14 were analyzed using all the data 
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points provl4«d In the 80 to 120 km region. The 
reaulte are 

03 ( 8 ) ■= 3.8200 X 10* m^'/kg 
Cjda) - 3.2352 X 10* m*/kg 
Cj(14) - 5.2345 X 10* mVkg 

The data for soundings 2 and 3 were noisier 
than the rest and a value of C 3 could only be 
obtained after eliminating a number of the moat 
divergent points. The results were 

03 ( 2 ) - 5.852 X 10* m*/kg 

€ 3 ( 3 ) - 15.924 X 10* m*/kg 

The noisy data associated with these results Is 
likely the cause for the values being higher than 
for the other three. 

The data for sounding 12 was smooch yet a 
solution for C 3 could not be obtained. It should 
be noted, however, chat sounding 12 was also cause 
for concern by the Michigan group because of the 
anomoulous behavior which they felt might be due 
to a small leak In Che Inflated sphere. 

Results From 7 In. Sphere 

Additional overlapping falling' sphere data was 
obtained for an accelerometer Instrumented 7 In. 
sphere experiment over Wallops Island flown In 1961.. 
The data Is published In NASA-CR-29 by Peterson and 
McWatters^ consisting of results obtained from Che 
first testa of the accelerometer system. The 
method of analysis was the same as used on the 
66 cm data. The result obtained was 

C 3 (7 In) - 1.5583 x 10* m^'kg 

This result Is within a factor of two of that 
obtained for soundings 8 , 13, and 14. The factor 
of two difference may be due to the different 
surface properties of the 7 In. as compared to the 
66 cm. Also, the 7 In. sphere enters the 
transition regime at a lower altitude than does 
the 66 cm sphere due to Its smaller size. The 
different molecular composition at lower altitudes 
may cause a change In C 3 . These questions could 
be answered by analysis of more 7 In. trajectory 
data. 

Atmospheric Density Calculations and 
Discussion of Results 

After determining a value for C 3 , the drag 
data taken In the original experiment may be 
reanalyzed to obtain new values of density employ- 
ing equation (22). An Iterative technique was 
employed using the unsmoothed temperature and 
drag values given in Reference 3. A comparison 
was made of the density values given In reference 
3 with the density values obtained using 


C 3 - 4 x 10*m^/kg which represents an average of 
the results for soundings 8 , 13, and 14. A point 
by point comparison was made for each of the 
soundings 8 , 13, and 14 and the average was 
obtained. The results are presented In Figure 2 
which shows that the densities calculated using 
the methods described In this paper produce 
generally higher values than obtained In the 
original experiments. 

Figure 3, obtained from reference 2, shows 
the mean values of density obtained In all 13 of 
the original experiments as compared to the 1962 
I). S. Standard.* This figure shows that the 
original analysis resulted In a nearly lOZ lower 
density value above 100 km than given In trefer- 
ence 1. Application of the results of the current 
work shown In Figure 2 would cause Che density to 
be nearly equal or somewhat greater than the 
U. S. Standard above 100 km. Both methods of 
analysis give a higher density In the 90 to 100 km 
region than that given In the U. S. Standard, but, 
neither analysis Is correct In this region as 
discussed In the following: 

In the 90 to 100 km region, the Knudsen number 
Is of order .1 as is shown In Figure 2. Transition 
flow Is considered to be within the limits of 10 
to .1 which means that in the 90 to 100 km region 
the flow Is near continuum rather than near free 
molecule. For this reason, equation (19), which 
Is derived on the basis of near free molecule 
theory. Is likely not valid In this lower region. 
The results could be Improved by employing a more 
valid relationship. The method of analysis would 
remain the same however. 

The results obtained In the 90 to 100 km 
region In the original analysis, reference 3, are 
also not correct due to inaccurate values of Cd> 

The recent work on CD reported In reference 6 shows 
that the values of Cjj used In the original 
analysis were -about lOZ higher than those measured 
In the wind tunnel. Therefore, somewhat higher 
values of density would be obtained In the 90 to 
100 km region but not to the degree Indicated by 
the results shown In figure , 2 . 

Above 100 km, an average lOZ higher density 
Is found using the Cp values calculated from the 
transition flow analysis. This difference can 
partly be explained by the difference In the 
treatment of the gas surface Interaction. In the 
current analysis, a value of K ~ 0 was used idilla 
In reference 2 the treatment of the gas surface 
Interaction resulted In an additive term of the 
form 2\/^3S„ was used where Is the speed 
ratio molecules would have If they obtain the 
temperature of the sphere wall, Tw « 300°K. Values 
of K which would compare more to the assumption 
made In reference 3 were attempted but the results 
for C 3 obtained for higher values of K were not as 
consistent as the ones obtained with K *• 0. In 
some cases, no solution for C 3 could be found for 
K greater than zero. These results tend to 
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indicate Chat K • 0 in a proper choice but more 
concluaive evidence ia needed in order to fix this 
important parameter. 

Conclusions 

The method of analysis reported here has 
demonstrated potential for application in high 
altitude falling sphere experiments. The values 
of C 3 obtained represent one of Che first experi- 
mental measurements of this quantity under high 
altitude conditions. More accurate values of C 3 
are needed to remove this unknown in Che analysis 
of falling sphere experiments could Chen be 
designed to measure still ocher unknown quantities 
such as temperature. Dual falling sphere measure- 
ments of both temperature and density at high 
altitudes have been proposed and appear to be an 
experiment well worth performing. 
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Figure 1. Velocity-altitude history for falling 

sphere sounding 13 in which ascent and 
descent trajectory overlap. 


48 



Altitudi (km) 


If 

1 



OERftfmmeiptRtM) 



Departure (percent) 


Figure 3. Departure from standard atmosphere of 
of the mean and standard deviations of 
13 density measurements deduced In the 
original experiments at Kwajaleln 
Island (plot Is copied from reference 2). 


Figure 2, Average of ratio of density computed 
using equation (19) and density 
obtained In original experiment for 
soundings 8 , 13, and 14 with 
C 3 “ 4 X 10* m^/kg. 
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CHAPTER V 


A SPHERE DRAG BRIDGING RELATIONSHIP IN THE LOW 
MACH NUMBER TRANSITION REGIME 
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A Sphere Drag Bridging Relationship In the Low 
Mach Number Transition Regime 
by 

** 

Gerald R. Karr 

The University of Alabama in Huntsville 
Huntsville, Alabama 35807 

Introduction 

Transition flow is defined as the flow regime in which the mean 

free path of the gas molecules is of the same order as a typical dimension 

of the body under consideration.^ At the boundaries of the transition 

regime are the free molecule regime (mean free path much larger than the 

characteristic dimension) and the continuum regime (mean free path much 

less than the characteristic dimension). The drag coefficient of a sphere 

is known to change markedly in the transition flow regime between the 

limits of free molecule and continuum flow and is a strong function of 

Mach Number. The drag coefficient in the free molecule regime approaches 

infinity at zero Mach Number and approaches a value near 2 at hypersonic 

Mach numbers. At the continuum limit, on the other hand, the sphere drag 

coefficient has a more complex nature which is known to depend on the 

Reynold's Number and the turbulence or lack of turbulence in the flow. 

In this work, however, the high density boundary of the transition regime 

4 

will be assumed to be at Re ~ 10 for which the drag coefficient has 

CO — 

a value of near 0.4 at subsonic velocities, increases in the transonic 
regime to a value of near 1.0 and approaches 0.92 at hypersonic velocities. 
Research supported under National Aeronautics and Space Administration 
Contract NAS8-28248, Marshall Space Flight Center, Huntsville, Ala. 
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4 

The continuum limit of the transition regime was taken as Re = 10 for 

00 

two reasons: (1) the sphere drag data employed in this work all corre- 
sponds to Re < 10^ and (2) sphere drag variations which occur above 

4 

Re = 10 are more clearly correlated with continuum parameters (Re and 
00 

turbulence) rather than what is normally considered transition flow 
parameters (Kn and surface to gas temperature ratio) . 

Due to the lack of adequate theory of the aerodynamics in the 
transitional regime, analytic determination of the sphere drag coefficient 
is usually made through semi empirical relations which are based on near 
free molecule flow theory and experimental results. These formulas are 
called bridging relationships, a number of which are reviewed in refer- 
ences 2, 3. and 4. The accuracy of a bridging relationship may be de- 
termined by comparison with experiment and most formula have at least 
one free parameter in order to obtain a best fit with given data. As 
discussed in reference 4, it is found, however, that available bridging 
relationships are typically accurate only over a limited range in 
Mach Number and Knudsen Number. The purpose of this paper is to report 
on a bridging formula which, with three free parameters, was found to 
predict to about 67 „ accuracy the sphere drag results obtained by the 
ballistic range method by Bailey and Hiatt. ^ Although Bailey and Hiatt 
provide plots of curve fits of drag coefficient as a function of Reynold' s 
Number to a claimed accuracy of ± 2%, the analytic results obtained here 
are provided as a function of Knudsen Number and have the advantage of 
allowing for ready interpolation as a function of Mach Number and Knudsen 
Number. Finally, it should be noted that the results obtained here are 
applicable only to fitting the Bailey and Hiatt measurements which are 

somewhat unique for sphere drag data in that the sphere surface temperature 
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was equal to the gas temperature (T /T = 1) . The results then have 

W CO 

application, for example, to falling sphere data analysis where T /T 1 

W 00 

but would not be applicable to wind tunnel data where typically T^/T »1. 

The Bridging Relationship 

The bridging relationship used in this work is a modification of 
that developed by Matting^ and also given by Rott and Whittenburg.^ 

Using a first collision, two fluid flow approximation. Matting obtained 
the result which can be written as 


C 


D 




o ’ E/Kn 
e 


( 1 ) 


where is the drag coefficient, is the continuum drag coefficient, 
is the free molecule drag coefficient, Kn is the Knudseh Number 
(defined as Kn = mean free path/ sphere diameter), and E is the free 
parameter. Equation 1 is seen to provide the correct limits as Kn is 
allowed to vary. At the free molecule flow limit, Kn -» oo, which gives 

4 

. At the continuum limit which in this work is taken as Re = 10 , 

» »fm 

Kn -♦ 0, which gives 


The limits are approached asymptotically 


which is what is observed experimentally. As will be shown, however, 
Equation 1 is not found to accurately predict the variation in the 
low Mach Number transition regime for any value of the free parameter E. 
Equation 1 is found to predict a much steeper variation in than what 
is observed in recent experimental results of Bailey and Hiatt. ^ 

In order to correct the failure of Equation 1, a second free 
parameter was introduced which was found to improve the accuracy con- 
siderably. The new form of the bridging relation is given by 

.X 


c 


(%M ■ \ ) 


E/(Kn) 


( 2 ) 
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where x is the new free parameter Introduced here. By raising Kn to a 
power, the steepness of the variation of with respect to Kn may be 
controlled, thereby better fitting the experimental results. 

Method of Determination of Free Parameters 
The values of free parameters are determined from a best fit to 
experimental data. The experimental data used here is that reported by 
Bailey and Hiatt^ which are obtained by the ballistic range method for 
which T /T =1 and covers a range in Mach Numbers from 0.1 to 6.0 and 

W CO 

a range in Reynold's Numbers from 15 to 50,000. Due to a lack of coverage 

in the transition regime at the lowest Mach Numbers, the data used in 

this work is limited to between M sa .72 and M ?a 6,0. Also, since only 

4 5 

6 of the 356” data points in the range have 10 < Re < 10 , a value of 

4 

Re = 10 has been used as the continuum limit. This range in flow 

CO 

parameters is of particular interest to falling sphere data analysis and 
also has application to satellite reentry and sounding rocket trajectories. 

Bailey and Hiatt provide tables of the experimental results 
arranged in gtoups of approximately the same Mach Number. For example, 

32 measurements of sphere drag coefficients were made for 1.45 < M <1.65 

” CO “ 

and 36 measurements made for 2.8 < M < 3.2 (see Table I for complete 

00 — 

list). For each measurement the values of M , Re , and C_ are given from 

CO 00 D 

which a Knudsen Number can be derived using^ 



where JL is the mean free path, d is the sphere diameter, and y = 1.4. 

The continuum and free molecule drag coefficients are assumed to 

8 

be functions of Mach Number only. The free molecule drag coefficient 
is given by 
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(1+K) 


(4) 


where erf Is the error function and where S is the speed ratio given 
by S = V/ ^2 RT which can also be expressed In terms of Mach Number 
as S = M ^y/ll The quantity K In Equation 4 Is a factor of order unity 
dependent on the gas surface Interaction. 

Since the evaluation of Equation 4 Is complicated by the presence 
of the error function, a useful expansion of Equation 4 for low and high 
values of S was employed In the analysis. The results of the expansion 
are 




8 

15 


S - 


8 

210 


'] 


(1 + K) 


: (S > 1.25) = 2 + (1 + K) 

°FM L 4S^J 


(5) 


( 6 ) 


which are accurate to better than .17. with respect to Equation 4. 

The continuum values of drag coefficient were obtained also from 
Bailey and Hiatt using values of C versus M for Re = 10,000 which, for 

D CD 

-4 

the Mach numbers of interest here, correspond to Kn ^ 10 . The ex- 

pression employed Is given by 

Cjj (M > 1.0) = .92 + .166/M - .366/M^ (7) 

c 

which is found to give an accuracy of at least 5% with respect to the 
experimental results. For Mach Numbers below 1.0, the following values 
were used 

Cjj (.1 < M < 1.0) = .403 - .142 M^ + .459M^ (8) 

c 
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Using Che relations given above, the data In a given Mach Number 
group were used to find the best values of x and E In a least squares 
sense. The least squares equation was written and the partial derivatives 
with respect to x and E were found. A computer program was developed to 
find the values of x and E which made the derivatives zero and thereby 
made the error a minimum. 

Results and Discussion 

For each Mach Number set tested, a root-mean- square (rms) value 
was computed and used as a measure of the accuracy of the fit for that 
set of data. The value of K required for the free molecule drag co- 
efficient value was found to have Influence on the results obtained. 

The RMS value for a given Mach Number was found to be improved if the 
value of K was taken to be a small negative number. Therefore, K becomes 
a third free parameter of the fitting process. 

The least square results showed Chat x and K have a Mach Number 
dependence while E is nearly constant. The x and K values were found to 
be nearly linear in Mach Number and approximated by the following rela- 
tionships, 

X = .399 + .016 M (9) 

K = -.002438 - .01842697 M (10) 

and the average value of E was found 

E = .212 (11) 

A set of X, K, E values are thus obtained over the full range of Mach 
Numbers under consideration. Equations 9, 10, and 11 were employed in 
Equation 2 and the results compared Co the ARO drag data. The RMS values 
that resulted are given in Table I along with the number of data points 
and the values of x, K, and E. The average RMS value of the 16 sets of 
data tested in Table I is .059. 
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Figure 1 shows a plot of vs Kn for six of the seventeen sets 

of data employed In the analysis. . The bridging relationship Is shown 

as the solid lines which are calculated based on the midrange Mach 

number of a given set of data. The figure Illustrates the success of 

the bridging relationship' In the transition regime and shows that much 

of the 67o rms error can be traced to scatter In the data which Is nearly 

± 107„ at some values of Kn. One failure worth noting, however. Is the 

tendency of the bridging relationship to underestimate the value In 
-3 -2 

the 10 < Kn < 10 range by about 5% In some cases. This Is likely 

a slip flow Influence which has not been taken Into consideration In 
this work. 


Table I 

RMS Values of Curve Fit to ARO Data Using Equations 9, 
10, and 11 in Equation 2 


Midrange Mach Number of 


Mach Number 

Number 

Range 

Data Points 

X 

E 

K 


.72 

.70 

- 

.74 

7 

.4105 

.212 

-.016 

.013 

.81 

.79 

- 

.83 

9 

.4120 

.212 

-.017 

.076 

.915 

.88 

- 

.95 

20 

.4136 

.212 

-.019 

.065 

.965 

.96 

- 

.97 

10 

.4144 

.212 

-.020 

.123 

.989 

.98 

- 

.998 

8 

.4148 

.212 

-.021 

.082 

1.135 

1.08 

- 

1.19 

24 

.4172 

.212 

-.023 

.036 

1.25 

1.2 

- 

1.3 

15 

.4190 

.212 

-.025 

.055 

1.375 

1.3 

- 

1.45 

28 

.4210 

.212 

-.028 

.075 

1.55 

1.45 

- 

1.65 

32 

.4238 

.212 

-.031 

.050 

1.75 

1.65 

- 

1.85 

23 

.4270 

.212 

-.035 

.051 

2.05 

1.9 

- 

2.2 

30 

.4318 

.212 

-.040 

.060 

2.55 

2.4 

- 

2.7 

18 

.4398 

.212 

-.049 

.059 

3.00 

2.8 

- 

3.2 

36 

.4470 

.212 

-.058 

.054 

A. 00 

3.8 

- 

4.2 

28 

.4630 

.212 

-.076 

.060 

5.00 

4.8 

- 

5.2 

40 

.4790 

.212 

-.095 

.041 

6.00 

5.8 

- 

6.2 

33 

.4950 

.212 

-.113 

.042 
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Conclusions 


A sphere drag bridging relationship has been developed for the 
low Mach Number transition flow regime which fits recent experimental 
results to an accuracy of about 6?, rms. The experimental results used 
were exclusively those reported by Bailey and Hiatt in tests ran at ARO 
and reported in March 1971 in which the sphere surface temperature was 
equal to the gas temperature. The results of this work should have 
application, for example, to the analysis of falling sphere data in 
which T /T « 1. 

W CO 

Due to the unique nature of the Bailey and Hiatt data (i.e., 

T /T = 1), it is of interest to examine the conclusions these data 

W oo 

indicate concerning the nature of the transition flow regime. Using 
the parameters found in fitting these data, Eq, 2 was plotted, vs M, 
(Fig. 2) for 0 < Kn < 00 and .1 < M < 6 for constant values of Kn. Since 
the highest Kn tested by Bailey and Hiatt was 10 ^ and since there was 
little transition flow data below M = 1, the curves for Kn > 10 ^ and all 
the curves for values of Kn for M < 1 are extrapolations. The results 
obtained, however, point to two important conclusions: (1) the width of 
the transition regime in terms of Knudsen Number is wider than usually 
assumed and (2) the free molecule drag coefficient implied by the results 
is less than usually assumed. 

Since the value of x (which is a measure of the slope of the 

vs. Kn curve in the transition regime) was found to be about .45, the 

width of the transition regime is increased over that obtained using the 

Matting relation which has x = 1 (Equation 1) . This conclusion is 

Illustrated by substituting into Equation 2 the value Kn = 5 which is 

the usually assumed upper limit of the transition regime.^ The results 
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at M = 6, Kn = 5 give C /C = .956 which shows that the free molecule 

^ °fm 

limit has not been reached at this value of Knudsen Number, A value 

of = .99 is found to be reached for M = 6 at a value of Kn = 

fm 

100 . 

The free molecule limit was found to be of importance in the de- 
velopment of the bridging relationship since K had to be adjusted to neg- 
ative values in order to obtain an accurate fit. This implies a free 
molecule drag coefficient less than two for the high values of Mach 

9 

Number tested. Results from other experimenters show that the drag 
coefficient in the free molecule limit is greater than 2 at Mach Numbers 
from 4 to 6. This departure from the results of others is likely ex- 
plained in that much of the sphere drag data at high Knudsen Numbers and 
high Mach Numbers used by others are obtained in low density windtunnels 
whereas the experimental data employed in this study was obtained ex- 
clusively from ballistic range data. The higher surface- to- gas- tempera- 
ture ratios that occur in windtunnels cause higher free molecule drag 

due to the energetic reflection of molecules at the surface. This effect 

9 10 

is clearly shown in recent sphere drag experiments ’ in which C 

^fm 

is found to decrease as T /T is decreased. In fact, an extrapolation 

W CO 

of the data of reference 10 down to T /T =1 indicates a C of 2 or less 

W oo J) 

A free molecule sphere drag coefficient of less than 2 requires that the 
gas surface interaction be non uniform on the surface of the sphere. 

This possibility is discussed by Cook^^ in which he proposes that the 
lowest possible free molecule drag of a sphere is 1,5. The lowest value 
implied in this work is 1.844 at M = 6 where K = -.113. 

The above conclusions are tentative but point to a need for free 
molecule experimental data in the low Mach Number regime at surface 
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temperatures close to ambient. Although the Bailey and Hiatt data do 
not reach free molecule conditions, the extrapolation discussed above 
show that the data lead to different conclusions than obtained from 
T /T ^1 experiments. 

W CO 
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FIGURE CAPTIONS 


Figure 1. Comparison with Equation 2 of Sphere Drag Coefficient Data 
(Ref. 5) as a Function of Free Stream Knudsen Number for 
Constant Values of Mach Number. 

Figure 2. Sphere Drag Coefficient as a Function of Free Stream Mach 
Number for Constant Values of Knudsen Number. 


62 




Figure 2 
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Improvements in Falling Sphere Data Analysis in the 80 to 120 Km Region 

by 

Gerald R. Karr and Robert E. Smith 

The analysis o£ falling sphere drag data is a principle means of 
density and temperature determination in the 80 to 120 Km region of the 
earth's atmosphere. This important method of atmospheric probing was 
reported by Bartman, et al^ in 1956 and has found wide use in upper atmo- 
spheric research. The method is particularly useful in the 80 to 120 Km 
region which is above the altitude capability of most aircraft and balloon 
probes but below normal satellite altitudes. 

In a typical falling sphere experiment, a sphere is ejected from a 
sounding rocket and the trajectory of the sphere is measured as it falls 
through the atmosphere. The trajectory information is analyzed to de- 
termine the velocity and acceleration as a function of altitude. This 
information is then used in the drag equation 

Drag = p Cjj A (1) 

where density, p , is obtained for a given value of drag coefficient, C^, 
and sphere crosssection area A. 

Temperature values are obtained from the density measurement using 
the hydrostatic equation 

dp = - g p d z (2) 

where g is the acceleration of gravity, z is the altitude, and p is the 
pressure. Equation 2 may be integrated between any two limits in altitude 
to give the pressure at z 
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p(z) = - 


g p dz + p (z^) 


( 3 ) 


f 


The temperature at z is obtained by substituting the values of p and p 
into the equation of state. 


p(z) = p (z) (R/W(z)) T(z) (4) 

where R is the universal gas constant, W is the molecular weight and T 
is the temperature. 

The temperature determination is seen to require knowledge of the 

pressure or equivalent temperature at some reference altitude (p(z ) in 

o 

Equation 3). In practice, the integration of Equation 3 is taken in the 
negative direction, from the point of highest ascent down to the lower 
altitudes. The temperature at is usually obtained from some atmospheric 
model. The error caused by possible incorrect temperature selection is 
minimized by the practice of downward integration since the term P(^q) 
in Equation 3 becomes small in comparison to the first term as the inte- 
gration proceeds to the lower altitudes. Thus, the effect of error in the 
initial temperature selection should become unimportant at one to two scale 
heights below the initial selection point. This was verified in the present 
study by selecting temperatures from 0°K to 1000°K at 120 Km. The effect 
on results below 100 Km due to such wide choices in temperature at 120 Km 
was found to be negligible ( < 1%). 

Another source of error in the falling sphere method is the value of 
Cjj used in the drag equation. Equation 1. The typical trajectory of falling 
spheres in the 80 to 120 Km region is found to correspond to an aerodynamic 

flight regime for which sphere drag coefficients have been highly uncertain. 
^ “5 *8 3 

Due to the low densities (10 to 10 Kg/m ) and the low Mach numbers 
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(1 to 5) experienced by falling spheres in this 'region, the aerodynamic 
flight regime Is classified as the transition flow regime (Ref. 2). The 
mean free paths of gas molecules in this region of the atmosphere vary 

_3 

between 10 m to lOm which means that the flow is too dense to be considered 
free molecule but two rarefied to be considered continuum. While the free 
molecule and continuum values of sphere drag coefficients have long been 
well known, the values of in the transition flow region have only re- 
cently been measured to adequate accuracy. 

The purpose of this paper is to report on improvements made in the 
falling sphere method of analysif. The most important of which is an im- 
proved relationship for the sphere drag coefficient which has estimated 
accuracy of at least 5 % over the range applicable to falling sphere tra- 
jectories. The second improvement in the analysis is to employ both the 
ascent and descent trajectory data in the data reduction. These improve- 
ments are applied to published falling sphere data and comparison is made 
with the results of previous analysis. 

Drag Coefficient Relationship 

3 

Recent sphere drag experiments reported by Bailey and Hiatt have 

provided and improved the drag coefficient information in the transition 

regime. The experiments were made in a Mach number and Reynold's number 

range of particular application to falling sphere data analysis (0.1 < 

M <6.2 and 20 < R < 10^). These data were obtained in ballistic range 
00 e 

in which the temperature of the gas and the sphere temperatures were 
approximately the same. The data are obtained in a consistent manner by 
one group of experimenters over the complete range of flow parameters 
applicable to the falling sphere flight regime. In view of the importance 
of this single source of drag data, a curve fitting relationship employing 
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this data exclusively has been prepared by Karr. The sphere drag bridging 
relationship thus developed is given. by 


c 


-E/a (Kn)* 




(5) 


g 

where a = .499 and where C and C are the continuum and free 

°c ®FM 

molecule drag coefficient, respectively. The quantities E and x are 

parameters of the curve fit and Kn is the free stream Knudsen number. 

The free molecule and continuum drag coefficients may be written as 

4 

functions of speed ratio, S , and Mach number, M, given by 

,3 




( 6 ) 


C (S > 1.25) = 2 11 + 

^FM 


-2 - 

V S 45 / 


(1 + K) 


(7) 


Cp ( M > 1.0) = .92 + .166 /m - .366/M^ (8) 

c 

Cjj (.1 < M < 1.0) = .403 - .142 m^ + .459M^ (9) 

c 

where 

S = V/ •v/2 RT ' and M = V/ ^ y El' (10) 

' W ^ W 

where y is the ratio of specific heats of the gas, and (1 + K) in 
Equations 6 and 7 is a factor or order unity which is related to the gas- 
surface interaction. The least squares curve fit of the above relations 
to the data of Bailey and Hiatt has revealed that Equation 5 will fit the 
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data to an accuracy of about 5% (based on RMS value) for the following 
values of the parameters 


X = 

.399 + .016 M 

(11) 

K = 

-.002438 - .01842697 M 

(12) 

E = 

.212 

(13) 


The use of Equations 6 through 13 in the sphere drag bridging rela- 
tion given in Equation 5 provides values of drag coefficients which are 
based on the most accurate and applicable Information now available. In 
order to use this relationship in the drag equation, values must be 
available for the velocity V, temperature T, molecular weight W, and the 
Knudsen number, Kn. The velocity is a measured quantity while values for 
the others are obtained as discussed in the following. 

Knudsen Number 

Knudsen number is defined as the ratio of near free path, to the 
characteristic length of the object. We take the characteristic length 
to be the sphere diameter, d . Therefore, 

Kn = i / d (14) 

The mean free path is inversly proportional to number density, n, for a 
simple gas (c.f. Ref. 5) and given by 

2 

i. = ( 2 TT n a ) (15) 

where a is the diameter of the gas molecules. Using the average molecular 
weight, we write the mean free path as a function of density 

‘ - (2 " m7^ 
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where is the Avogadro constant. Using a value of a = 3.72 x 

which is representative of air^, we obtain 

Kn = 2.7 X 10’^ (17) 

3 

for p in units of Kg/m and d in meters. 

In developing Equation 17, we have neglected the small effect tempera- 
ture has on the value of a and we have assumed the gas molecules each 
have a mass corresponding to the average molecular weight of the gas. 

Molecular Weight 

A value of molecular weight is required for the Mach number, speed 
ratio, and Knudsen number expressions. (We assume y ~ independent 

of the molecular weight) . The mean molecular weight is known to be a con- 
stant of 28.964 up to about 90 Km altitude. At that altitude, the heavier 
molecules begin to settle out causing a drop in molecular weight. Models 
of this molecular weight with altitude are provided in Ref. 6. Representa- 
tive values of the variation is given by the following equation 

W = 24.68 + ,1235 z - .000874 z^ (18) 

where z is the altitude. Equation 18 gives values of molecular weight 
with an accuracy of better than 17o in comparison with a nominal spring/ 
fall values given in Ref. 6, 

Temperature and Density Iteration 

The temperature at a given altitude is determined as stated in the 
introduction, by downward integration of the density profile. However, 
since the drag coefficient is a function of temperature (i.e., Mach number 
and speed ratio are inversely proportional to the square root of temperature), 
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values of temperature xnust be available before determination of density 
is made using the drag equation. Thus, we are lead to an iterative 
procedure for finding both density and temperature. The procedure is 
begun by assuming some Initial temperature profile from which, given the 
unusual velocities, the Mach number and speed ratio are found. Once a 
density profile is obtained (to be discussed in the next paragraph) using 
the assumed temperature profile, a new temperature profile may be con- 
structed. This process is repeated until the values of both temperature 
and density no longer change beyond a specified error limit. 

Since the drag coefficient has been written as a function of density, 
the drag equation becomes of the following form 

X 

Drag = ~ p A C + e (C - C ) (19) 

c FM - ■ 'c 

where b is defined through Equation 17 and 5. 

b = E/(2.7 X 10"^ W/d)^ (20) 

Equation 19 is a nonlinear equation in the unknown density and must be 
solved using numerical procedures. 

Solving for density from the drag equation of the form given in Equation 
19 has two advantages. First, as already pointed out, the relationship 
for Cjj represents recent, accurate drag measurements. Second, by solving 
Equation 19 for density directly, we eliminate the uncertainty of choosing 
a drag coefficient from a set of values tabulated as a function of Mach 
number and Reynold's number. Falling sphere experimenters effectively 
solve an equation like Equation 19 by iterating between tabulated values 

and using the simple drag equation. The solution of Equation 19 is obtained 
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much faster and the solution is likely more accurate than that obtained 
in such iterative procedures. 

Application of the Proposed Method 

The procedure outlined above was applied to sets of falling sphere 
data reported by Peterson, et al^ for measurements made over Kwajaleln 
during 1963 and 1964. One of the reasons for using this data to illustrate 
the application of the proposed method of analysis is the importance of 
the results which were obtained in the original experiment. This set of 
data is of particular interest since it forms the basis for the density 
and temperature model in the 80 to 120 Km region for what is labeled 
"15°N, Annual" in the U. S. Standard Atmosphere Supplements, 1966.^ A 
summary of the density results obtained by Peterson, et alj are given in 
Figure 1 which is a copy of Figure 2.20 of Reference 6. 

A second reason for using this particular set of data was the avail- 
ability of the information required in our analysis. Through one of the 

g 

original experimenters, K. McWatters, we obtained the detailed computer 

output (a sample of which for one flight is given in Reference 7) for the 

falling sphere measurements made by the group. The data of interest is 

the density, p^, temperature T^, drag coefficient Reynold's number 

Mach number M , and velocity V , as a function of altitude, z , which 
o’ ft’ n 

resulted from their analysis of the falling sphere trajectory. Of these 
quantities, only velocity and altitude were measured where the subscript 
n is used to indicate this. The remaining quantities were deduced and we 
use the subscript o to indicate this. The measured drag force, however, 
can be found by taking the deduced density and the drag coefficient em- 
ployed in obtaining that density, and substitute these values into the drag 
equation. Therefore, 
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(Drag)^ 



where the area is based on 0.66 m sphere diameter. Thus, in this manner, 
the values of altitude, velocity, and drag force are obtained as measured 
in the original experiments. The accuracy of these measurements is ex- 
pected to be at least 5%? 

In order td begin the iteration for density and temperature, the data 

analysis requires that an initial temperature profile be given. Of all the 

temperatures given in this initial profile, the only one of importance is 

the one given at the altitude at which the downward Integration is begun. 

For the Kwajalein data, the highest altitude for which data are available 

is 120 Km since data above 120 was considered too inaccurate. In our 

analysis we have taken the temperature at 120 Km to be 460°K which is based 

on the 15°N Annual model given in Reference 6. For comparison, Peterson, 

et al? used a temperature value of 361°K at 120 Km which is based on the 

9 

1962 U. S. Standard Atmosphere. The higher temperature was chosen in this 
work in view of the fact that the results obtained in the original analysis 
revealed that the temperatures were generally higher than the USSA 1962 at 
the high altitudes. The effect the choice of a 100°K higher temperature has 
on the final results will be discussed more in the conclusions. 


Temperature Interpolation for the Descent Trajectory 

At this point in the discussion, we describe another improvement to 
the data analysis which is applicable to the Kwajalein falling sphere data. 
Much of the data obtained over Kwajalein consists of two trajectories; 
ascent and descent. Soon after ejection of the sphere, the radar tracks 
the sphere during its ascent to apogee. This tracking begins usually near 
100 Km and data is then terminated at 120 Km for accuracy reasons as dis- 
cussed above. During the passage through apogee, the radars "lose" the 
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sphere. Also because of smoothness requirements, data processing during 

(See Fig. 2) 

descent generally does not start again until about 100 Km. In the orig- 
inal analysis of this data by Peterson, et al7, the ascent and descent 
trajectories were analyzed separately. That is, the downward Integration 
required for temperature determination was begun again at the top of the 
descent trajectory with a temperature choice based on USSA 1962. The 
temperature thus chosen at the top of the descent trajectory is potentially 
erroneous, the effect of which will propagate two scale height down into 
the descent trajectory. We have developed a more accurate procedure for 
choosing the temperature value at the top of the descent trajectory which 
employs the ascent temperature values and an isentropic relationship. The 
ascent temperature values in the 100 Km region are considered to be more 
accurate than the higher altitude values since this corresponds to two 
scale heights below the 120 Km altitude where the temperature is taken 
arbitrarily to be 460°K. 

The temperature determination procedure employed in the present 
analysis is as follows; 

Step 1. Density profiles are determined for both the ascent and de- 
scent trajectories, p. (z) and p (z), respectively. The 15° N Annual model 

A U 

is used to give an initial temperature profile. 

Step 2. A new ascent temperature profile, T (z), is determined for 

A 

the ascent data using the downward integration method. 

Step 3. For the region of altitude for which both ascent and descent 
data points are available, an isentropic relationship is employed to find 
the temperature for the descent trajectory points, Tp(z), given by 

Y - 1 

Id(») - (pb(') / Pa<“>) 
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where y Is the ratio of specific heats. In using this relationship, we 
are assuming that the densities which have been measured at the same 
altitude by the drag method are different due to an Isentroplc process. 

The sphere trajectory Is such that at about 100 Km altitude, the ascent 
and descent trajectories are about 80 Km apart. Wave motion In this 
region of the atmosphere could then account for the differences In density 
which are seen. Over the length and time scales of Interest here, the 
assumption that such processes are Isentroplc Is reasonable. 

Step 4. For the remaining altitude points In the descent tra- 
jectory, the temperatures are determined using the downward Integration 
method. 

Step 5. Based on the new temperature values, calculate new 
values of Mach number and drag coefficient to be employed In the 
determination of new density profiles. That Is, start again at Step 
1 above and continue the Iteration procedure until convergence Is 
reached In both temperature and density results. 

The above procedure requires at Step 3 that at least one data 
point of the descent trajectory be at the same, or nearly the same, 
altitude as a point on the ascent trajectory. This requirement Is met 
for six of the thirteen flights made over Kwajaleln and reported In 
reference 7. Table I gives some of the characteristics of the flights 
studied In this work Including the range of overlap for the six flights. 
Two of the six flights had only one point In common while flight # 13 
was found to have 12 data points In common (data Is provided approximately 
every kilometer, on the kilometer). 
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Results 


Density and temperature profiles were obtained by the methods de- 
scribed above for the six falling sphere flights over Kwajaleln which 
had overlap In altitude coverage of the ascent and descent trajectories. 

The density profiles were compared with that published In the U. S. 

9 

Standard Atmosphere 1962. The ratio of the density found In this work 
to the USSA 1962 densities Is shown for each of the six flights In Figures 
3 through 8 . Also contained In these figures are the temperature pro- 
files obtained In the analysis. Table 1 lists the sounding number which 
was designated by the experimenters, the time and date of the flight, 
and the altitude range covered on the ascent and descent trajectories for 
the six flights analyzed. 

The mean and standard deviation of the density and temperature pro- 
file ratios (with respect to the models indicated) are plotted in Figures 
9 and 10 , respectively. For the altitudes outside the region of over- 
lap, six values were used in construction of the mean and standard de- 
viation, For some altitudes within the overlap region, as many as 12 values 
were used. 
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Discussion of Results 


The summary of density results given In Figure 9 show a signifi- 
cant departure from that obtained In the original analysis of Peterson, 
et al. 1966. This departure Is found to be primarily because of the In- 
complete Information on drag coefficients available to the original 
experimenter at that time. The drag coefficient values employed In 
the present analysis are considerably more accurate and reveal a point 
of maximum departure from the 1962 standard at 105 Km rather than at 
92 Km. 

The temperature results given In Figure 10 show better agreement 
with the original analysis than does density. The departure at 120 Km 
Is artificial since the choice of a temperature value at 120 Km Is 
completely arbitrary. The departure in temperature between the present 
and original analysis at high altitudes Is not the cause for the de- 
parture in density at these altitudes. In fact, if the value of 
temperature at 120 Km used by Peterson, et al. were to be employed in 
the present analysis, the departure in density results would be even 
greater than given in Figure 9. The colder temperatures used by 
Peterson, et al. would result in lower drag coefficients and therefore 
even higher densities would be obtained from the drag equation. The 
results we have obtained seem to clearly point towards the higher 
values of temperature at altitudes near 120 Km. 
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Table 1 


Figure Number 

Sounding Number 

Time 

(GMT) 

& Date of Flight 

Range of 
Ascent 
Altitude 

Range of 
Descent 
Altitude 

Overla 

Range 

3 

2 

0257 

March 29, 1963 

100-120 

102-80 

100-10 

4 

3 

0328 

June 18, 1963 

99-120 

102-80 

99-10 

5 

8 

1458 

Nov. 14, 1963 

104-120 

109-80 

104-10 

6 

12 

1820 

March 13, 1964 

99-120 

104-80 

99-10 

7 

13 

1125 

May 12, 1964 

96-120 

107-80 

96-10 

8 


0101 

June 17, 1964 

102-120 

102-80 

120 




HP 





Departure (percent) 


Fig. 1. Departure from standard atmosphere of the mean and standard 

deviations of 13 density measurements deduced in the original 
experiments at Kwajalein Island (plot is copied from Ref. 2). 



Fig. 2. Velocity-altitude history for falling sphere sounding 13 
in which ascent and descent trajectory overlap. 
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Temperature, K Density Ratio, p/p (ussa,i 962 ) 


Figure 3. Temperature and Density Profiles Obtained fron 
Sounding # 2. 



Temperature, K Density Ratio, /=/^( ussa,i 962 ) 


Figure 4. Temperature and Density Profiles Obtained from 
Sounding # 3. 
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Altitude, km Altitude, km 


II 



Temperature,K Density Rdtk),/?<^(ussA,i 962 ) 


Figure 5. Temperature and Density Profiles Obtained from 
Sounding # 8. 



Temperature , K Density Ratio, /o//> (ussA.isea) 


Figure 6. Temperature and Density Profiles Obtained from 
Sounding # 12. 
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Altitude, km 



Temperature, K Density Ratio, P/P (ussa.i 962 ) 


Figure 7. Temperature and Density Profiles Obtained from 
Sounding // 13. 



Temperature , K Density ^X\o,p/p (ussa,i 962 ) 


Figure 8. Temperature and Density Profiles Obtained from 
Sounding // 14. 
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ALTITUDE, Km 



DENSITY RATIO, ///“CUSSA, 1963) 


Figure 9. Mean and standard deviation of density results obtained 
in the present analysis for six falling sphere flights 
over Kwajalein. Mean of results from original analysis 
given for comparison. 


83 




TEMPERATURE/TEMPERATURE(USSA,©623 


Figure 10. Mean and standard deviation of temperature results 
obtained in the present analysis for six falling 
sphere flights over Kwajalein. Mean of results from 
original analysis given for comparison. 
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CHAPTER VII 


FREE MOLECULE DRAG AT SPEED RATIO LESS THAN UNITY 



Free Molecule Drag at Speed Ratio Less Than Unity 

by 

Gerald R. Karr 

Assistant Research Professor 
Mechanical Engineering Department 
The University of Alabama in Huntsville 

This work concerns results obtained in the calculation of the 
drag force acting on objects which (1) have dimensions much less than 
the mean free path of the gas (Knudsen number much less than one) and 
(2) have velocity with respect to the gas which is much less than the 
thermal velocity of the gas (Speed ratio, S, much less than one). These 
conditions are characteristic of those experienced by certain aerosols 
and Brownian type particles which are suspended in a gas or which are 
forced to travel through a gas such as in separation processes. 

The purpose of this work is to investigate the influence of two 

factors which enter the calculation of the free molecule drag force; 

(1) the shape of the object and (2) the gas surface interaction. Past 
12 3 

investigations ’ ’ have concentrated on the second of these factors 
while assuming the particles are perfectly spherical in shape. The 
gas surface interaction has generally been taken to be composed of a 
specular fraction and a second fraction which is purely diffusive (see, 
for example. Ref. 4) or a modified diffusive such as the elastic- 
diffusive reflection employed in Ref. 3. The form of the diffusive 
fraction of the reflected molecules has received considerable attention 
because comparison of sphere drag calculations with experiment lead to 
the conclusion that all the molecules must be diffusively reflected in 

order to explain the observed high drag coefficients at low speed ratios. 
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In the work reported in this paper, the drag coefficients for both 
specular and diffuse reflection were obtained for non-spherical objects 
in order to investigate the influence of shape at very low speed ratios. 

The drag coefficients were obtained by employing the expression 
for force acting on an element of surface in a free molecule flow valid 
at any speed ratio and taking the component of that force which acts 
-in the direction of the velocity vector. The drag component of force 
was divided into a part due to the momentum of impinging molecules, 

D^, and a part produced by the reaction force of the molecules leaving 
the surface, D^. Expressions for both specular D^(specular) and 
diffusive D^(diffusive) reflections were developed. The total drag 
force is then given by D = D^. + for an element of surface at any 
angle with respect to the flow and any speed ratio. The expressions 
were then integrated over the surfaces of various shapes including 
oblate spheroids, cylinders with flat and spherical ends, and cones 
with flat ends. The drag force coefficients were obtained for these 
objects at low speed ratios for both specular and diffusive reflection. 
The total drag force acting on an object is written 

D = D. (1 + ^ ) 

1 

where the ratio D^/D^ represents the effect of the reflection normalized 
with respect to the incident contribution. For a perfect sphere, for 
example, and for S « 1 

(sphere) = p tt r^ ^2 tt RT' V 

The attached figure presents examples of values of for oblate 

spheroid shapes as a function of the minimum-to-maximum radius ratio. 
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The results show the sensitivity of the tb both the reflection 

characteristic and the object shape. 

The results obtained in this work show that the magnitude of 
the force coefficient is strongly dependent upon the shape and o*"ien- 
tation of the object for both specular and diffusive reflection. Since 
non-spherical aerosol or Brownian particles would present randon orien- 
tation with respect to the velocity vector, the force coefficient for 
a given object would be an average value. The results obtained in this 
work provide the information needed to obtain the average force co- 
efficient for various non-spherical shapes. One conclusion reached in 
this work is that both specular and diffuse reflections can produce 
high drag force coefficients at low speed ratios for non-spherical 
objects and neither should be excluded from consideration in such cases. 
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A 


Proposal to Develop Zero-g Brbwnian Motion Experiments 

by 

G. R. Karr 


Introduction 

Robert Brown in 1828 is credited with establishing as an impor- 
tant phenomenon the observed irregular and perpetual motion of macro- 
scopic particles suspended in gases or liquids. The theory of this 
motion, which has become known as Brownian motion, has received the 
attention of Einstein, van Smoluchowski, Langevin , Uhlenbeck and 
Ornstein, Chandrasekhar, and many others. The motion received much 
early interest because it has been established through the theory that 
Brownian motion is direct observation evidence of the molecular state 
of matter. For example, from observation of Brownian motion, one can 
measure Avrogodros number and this method was in fact considered to 
provide the most accurate measurement of this quantity in the early 
1900' s. The theory of Brownian motion establishes that the motion is 
described as a random process which is found to be a major step in the 
development of the field of study now called stochastic processes. 

Experimental investigations of Brownian motion has not received 
the interest of physicists recently and modern interest in Brownian 
motion is primarily in the theory and its application. However, with 
the unique environment which will be provided by the Space Shuttle, the 
possibility now exists for performing experiments involving Brownian 
motion that could not be done in a one-g environment. It is proposed 
here to study possible experiments that may be performed in a zero-g 
environment which involve the observation of Brownian motion. It is pro- 
posed to assess the feasibility and importance of such experiments. 
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The Theory of Brownian Motion 


The currently accepted theory of Brownian motion was presented 
by Ornsteln and Uhlenbeck In 1930. Beginning with the equation of 
motion of a Brownian particle (called the Langevln equation) we have 

d^x 

m — r = - m p V + F(t) (1) 

dt 

where g Is a drag coefficient and where F(t) Is a random forcing term 
characteristic of the Brownian motion. Ornsteln and Uhlenbeck obtain 
the solution for the mean square displacement of the particle given by 

^ Ot - 1 + e ■ P*") (2) 

ra 3 

which has the limits for t large compared to p ^ 

which Is the result obtained by Einstein and for t small compared to 
P ^ the result Is 



(4) 


where u^ Is the Initial velocity of the particle. 

Other quantities may also be calculated. For example, the mean 
square velocity of Brownian particles all starting at velocity u^ is 




2 ^ t 


(5) 


and the velocity distribution function of the particles Is given by 
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GCUq, 


u. t) 


m 


TT k T(l-e 


-2 3 


( 6 ) 


which shows that the particles eventually reach a Maxwell-Boltziiiann 
distribution. . . - 

The above results are all for what is called a free particle. 

That is, the only forces acting on the Brownian particle are the. drag 
term m3v and the random forcing term F(t). For the case of Brownian 
motion in an external force field, such as gravity, theory is not so 
clear nor as complete as for the free Brownian motion. Uhlenbeck and 
Ornstein, for example, consider the Brownian motion of a particle which 
is bound in a harmonic force field at frequency Three cases result 

and the solutions for the mean square displacement under these conditions 
are provided. 

Overdamped case; 6 » 2 u) 


nu) ' m i JU / 


8 2 

where iju' = — - uj 

Critically damped case; S = 2 uj 


6 ' 
(cosh lo't + — ^ sinh uj't)' 

2 U)' 


(7) 


2 o 


kT 


m ju 

Underdamped case; B < 2 U) 

kT . / 2 kT 




2 ^o 
X m 


T ^ / 2 kT \ -Pt / 

“2 + 1^0 2 ) ® '"l 

d) \ m uj / \ 


t + sin U) 

2u)i 




U 2 2 

where = uj 


4 


( 8 ) 


(9) 


While Uhlenbeck and Ornstein obtained solutions for the harmonic 
forcing cases, the case for constant forcing has not been solved as 
completely. 
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For example, Wang and Uhlenbeck state In their 1945 paper that they 
have been unable to find the general solution for the constant force 
case. Solutions do exist, however, for the special case when the 
friction force is high and the observation time Is large (l.e., for 
t » p ^). Chandrasekhar, for example, obtains the time varying dis- 
tribution of Brownian particles in a gravity field and is able to show 
that the particles arrange themselves in a barometric distribution with 
respect to the gravity vector. 

Solutions for other observables of Brownian motion in a constant 
force field are apparently not available. For example, the mean square 
displacement or the velocity distribution function in a gravity field 
were not found in the references cited. Part of this study will be di- 
rected toward a thorough investigation of the recent literature on this 
subject. From the literature that has been searched, however, it is 
apparent that Brownian motion in a gravity field will be much more com- 
plex than the free motion and that sensitive experiments in one-g would 
be strongly influenced by gravity. 

Zero-g Experiments 

The zero-g environment offers at least three advantages in the 
performance of Brownian motion studies. 

1. The theory of Brownian motion in zero-g is well established 
while the motion under gravity is difficult at best to interpret. 

2. Convection currents can be minimized or eliminated in zero-g. 

3. In zero-g, particles of larger sizes and masses can be em- 
ployed in Brownian motion studies. 

The purpose of the proposed work is to study possible Brownian 

motion experiments and to assess the value and feasibility of such 

experiments. While it is expected that other experiments will be proposed 
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and studied during this work, the following set of experiments have 
already received preliminary consideration and will be described 
briefly. 

A. Particulate drag In the transition regime 

As Is well known. Brownian motion Is observed In both liquids 
and gases. The drag or friction coefficient, p. In these media Is, 
however, considerably different. The friction factor In a liquid for 
example Is In the fluid dynamic flow regime called "Stokes f lot/' and 
Is described as a hlghly-vlscous flow. For a Brownian motion In a 
gas, however, the flow Is free molecule If the mean free path of the 
molecules Is large compared with the size of the particle. It Is 
therefore proposed that a possible zero-g experiment Is to vary the 
properties of the suspension medium or the particles, so as to obtain 
friction factor Information In the transition flow regime between the 
Stokes regime and the free molecule regime. The proposed work would 
consist of determining the range of particle and medium properties 
that would be required to probe the transition regime. Also to be studied 
Is the limits Imposed on such an experiment by the one-g environment. 

The transition flow regime has proved very difficult to probe In ground- 
based experiments and becomes progressively more difficult as the speed 
Is reduced. The experiments proposed here would provide data at the very 
low end of the velocity spectrum, a region for which little or no data 
now exists. This possibility of Increasing the range of understanding 
of fluid dynamic drag appears to be of great value. , 

B. Gas surface Interaction 

Brownian motion observations provide an excellent basis for 

studying the Interaction of molecules with surfaces. The motion Is a 

direct consequence of the bombardment of the particle surface with 
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molecules of the surrounding medium. The free molecule regime Is best 
for such observations since the motion Is Isolated from the effects of 
viscosity that occur In the Stokes flow regime. The proposed experi- 
ments would consist of ( 1 ) the. preparation of particles of known com- 
position and surface properties, ( 2 ) the preparations of gases of known 
composition and temperature, and (3)' collection of data on the resultant 
mean displacement and velocity distribution . to determine the friction 
factor 3 . The value of 3 determined in such an experiment can be di- 
rectly related to the effect of the gas surface Interaction. The data 
can then be correlated with respect to the gas and surface properties. 

The physics of the gas surface Interaction Is not well understood 
at present and experimental data of the type proposed here would be valu- 
able in identifying the important characteristics of the interaction. 

A factor of two variations in the friction coefficient, p, is theoreti- 
cally possible due to the effect of the gas surface Interaction. Ground 
based experiments such as the oil drop experiment do not have the sensitiv- 
ity needed to determine these effects. 

C. Verification of Brownian Motion Theory 

The Brownian motion theory proposed by Uhlenbeck and Ornstein is 
the currently accepted theory of the motion. The theory of Einstein and 
Smoluchowskl Is found to be the limiting case of the Uhlenbeck and 
Ornstein theory for large times (t » p for the free particle case. 
Edward Nelson further proposes that the Elnsteln-Smoluchowskl theory Is 
also limited to the case of large friction (p large) and points to ex- 
perimental results of Kappler and of Barnes and Silverman to show that, 
for the harmonic forcing case, the Elnsteln-Smoluchowskl approximation 
is invalid for the underdamped case even for t » p For the same 

reasons as mentioned above, the presence of the one-g field requires 
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that particles be small and light. For such cases then, it is difficult 
to probe the Brownian motion for short times (t < B"^). The motion for 
^short times can be used to verify .the Uhlenbeck-Ornstein theory. The 
results for t » p ^ are the same for both the Uhlenbeck-Ornstein and 
the Elnsfeln-Smpluchowski theories and thus such results can not be used 
to distinguish the theories. Due to the difficulty in making measure- 
ments in the short time on earth, there appears to be little, if any, 
experimental verification of the Uhlenbeck-Ornstein theory. It is, 
therefore, proposed that since the zero-g environment offers the 
opportunity to adjust the size of p over a wide range, a properly de- 
signed Brownian motion experiment in space would allow for the verifi- 
cation of theory of Brownian motion. The proposed work will seek to 
establish the conditions required to perform such an experiment. 

Proposed Work 

The three experiments proposed above are clearly of great value 
and preliminary study Indicates that they are also feasible. During the 
proposed study, these experiments and others will be considered in detail 
to determine the value of the experiment, the reasons that zero-g are 
required, and the feasibility of performing the experiment. Extensive 
literature searches and personal contacts will be made to ascertain the 
state-of-the-art knowledge of Brownian motion theory, motion experiments, 
friction coefficients, gas-surface Interaction, etc. Experimental tech- 
niques will be surveyed and recommendations made for space experiments. 

It is also expected that preliminary experimental development will be 
undertaken to test observation methods, data collection methods, and 
data analysis techniques. 
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B 


An Experiment Using the Molecular Beam Apparatus Proposed for Space Shuttle 
Title; Aerodynamic Force Measurements in Space 

Principal Investigator; Dr, Gerald R. Karr, PKe University of Alabama in 
Huntsville, Huntsville, Alabama. 

Summary of Proposed Work 

The feasibility is to be determined of employing the satellite orbit 
environment in the measurement of forces resulting from the interaction of 
atmospheric gas with solid surfaces at satellite velocities. To be con- 
sidered is an experiment designed to measure the aerodynamic forces acting 
on surfaces exposed to the high-velocity, low-density gas flow which is 
generated as a satellite travels through the upper atmosphere. In partic- 
ular, the use of the proposed Molecular Beam Laboratory will be considered 
for providing the required beam definition and orientation. Engineering 
and scientific gains would be generated by the results of this experiment 
which utilizes an aspect of the orbital flight environment not easily re- 
produced in ground based facilities. The study would evaluate means for 
measuring the aerodynamic forces on a selection of surfaces having a broad 
range of material and physical properties. The study would determine the 
desirable number of surface materials, the range of surface temperatures, 
the range in degree of surface contamination, the number of surface coatings, 
and the angles of attack to be tested in the proposed experiment. Finally, 
the feasibility would be determined of correlating the force measurements 
with changes in gas properties. 

Justification 

The determination of the feasibility of the proposed experiment is 

desirable .in view of the potential benefits the experiment would provide. 
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At present, satellite aerodynamic properties cannot be predicted accurately 
because of the lack of Information the proposed experiment could provide. 
Analytic studies reveal that satellite aerodynamic properties are a strong 
function of the character of the force caused by the gas surface inter- 
action. The gas surface Interaction, in turn, is expected to be a strong 
function of surface properties and surface orientation. The design of 
satellites to take advantage of (or to reduce) the aerodynamic forces and 
torques has not been possible because of the lack of Information on the 
forces caused by the gas surface interaction. Knowledge of the character 
of the surface forces and the major influences on these forces is necessary 
to the design of satellites to have specified aerodynamic drag, lift, and 
torque characteristics. Such knowledge is also needed in order to inter- 
pret the dynamic response of satellites In the atmosphere as is done in 
the determination of atmospheric density from satellite drag measurements. 

In addition to the engineering Information provided by the proposed 
experiment, the results would also contribute to the scientific knowledge 
of the gas surface interaction at satellite velocities. The expected 
results would compliment both orbital and ground-based molecular beam 
studies ^Ich provide force information indirectly. The proposed experi- 
ment would then serve to guide future investigations into the more subtle 
details of the Interaction. 

Finally, the study of the feasibility of the proposed experiment 
is justified on the basis that the experiment may be a relatively in- 
expensive method of obtaining valuable Information. Consequently, the 
experiment could require few equipment components with low development 
cost and short development time. The information gained would be of 
immediate engineering value and would be a valuable complement for future 
gas-surface experiments. 
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Method 


The objective of this study is to determine the feasibility of 
performing a space experiment to measure the aerodynamic forces oh sur- 
faces as a function of gas and surface properties and surface orientation. 
The study is divided into three areas: (a) Methods of making measurements, 
(b) Selection of surfaces and surface conditions, and (c) Correlation of 
results with gas properties. 

(a) Measurement Techniques 

The feasibility of making the measurements required will be 
investigated taking into consideration the expected low level of force 
and the perturbing Influences of environmental factors. To be considered 
is the feasibility of the aerodynamic forces acting on flat or shaped 
surface samples exposed to the gas flow generated by the motion of the 
satellite through the atmosphere. The perturbing influence of molecules 
reflected from the satellite may require that the surface samples be ex- 
tended on a boom ahead of the vehicle. Methods will be evaluated for 
measuring the forces, orienting the surface samples, and changing the 
surface properties. 

The accuracy of the proposed measurements is to be evaluated 
considering perturbing environmental Influences such as upper atmospheric 
wind and density fluctuations. Methods of calibration and monitoring of 
the environment will be considered as means of increasing the accuracy of 
the proposed experiment. 

(b) Selection of Surfaces 

The selection of surfaces to be tested in the proposed ex- 
periment will consider the need to reduce satellite payload weight and 
volume while yielding results of the widest possible interest. The se- 
lection of surfaces will be on the basis of providing Information of the 
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many factors thought to Influence the forces caused by the gas surface 
Interaction. Among the factors of Interest are the influence of surface 
material, surface temperature, surface roughness, surface coatings, sur- 
face contaminates, and surface angle-of-attack to the flow. The feasibility 
study would establish a series of experiments which best Isolate the in- 

i 

fluence of the individual factors,. The surfaces selected will span those 
used in satellite construction so as to provide engineering information 
for future design. 

(c) Measurement of the Influence of Gas Properties 

Since the gas surface interaction is influenced by both the 
surface properties and the gas properties, the feasibility is considered 
of determining the influence of the gas properties on the surface forces. 

The upper atmospheric gas composition, temperature, and degree of ioniza- 
tion are a strong function of altitude, geocentric latitude and longitude, 
and time. To be investigated is the possible correlation of the measured 
forces with the changes in gas properties that occur naturally over the 
orbit. The feasibility will be studied of identifying the gas-property 
influences on the forces caused by the gas surface interaction. Study 
will be made of the orbital parameters which provide the best conditions 
for the experiment. 

Personnel 

The principal investigator of the proposed study is Dr. G. R. Karr 
(resume attached) who has done considerable work on the theory of the gas 
surface interaction and satellite aerodynamics. Since Dr. Karr has pri- 
marily theoretical experience and capability, there is a recogtllzed naed for 
cooperation with personnel who have experimental capability and experience. 
In view of the good working relationship which exists between Dr, Karr, 

The University of Alabama in Huntsville, and NASA Marshall Space Fligiht 
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Center, Huntsville, it is proposed that the theoretical expertise of 
Dr. Karr be complemented with the experimental expertise of MSFC personnel 
such as Dr. P. Peters (a surface physicist) and/or Dr. R. Smith (an atmo- 
spheric physicist) both of the Space Science Laboratory at MSFC. 
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APPENDIX 


Listing of computer programs developed and employed under NASA Contract 
Number NAS8-28248. 


Program Name Page 


LESQA 

106 

RHORAT 

111 

AFILIP-HIGH CM 

119 

AFILIP-LOW CM 

123 

CDCLEV 

127 

CFEVAL 

130 

RUFSPH 

133 

CLELAN 

137 
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PROGRAM LESQA 


This Is a computer program which analyzes density, temperature, 
velocity and altitude measurements from falling sphere experiments. 
Also In the Input data are the densities and drag coefficients em- 
ployed by the original experimenters so that acceleration data can 
also be deduced. Both ascent and descent data are employed. The data 
are employed in an orthogonal curve fitting routine so that ascent and 
descent data can be correlated at equivalent altitudes. Speed ratio 
effects, molecular weight changes and drag coefficient variations are 
all taken Into consideration. A density Is determined which best 
represents the data based on the measured properties and those pre- 
dicted by theory. A density profile is thus determined using both 
ascent and descent Information. 
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I 


«RUN»/lP LtSQA»UAHXXXXXXXXX»tNVIH-DYN-faK»t5»lbO 


XXMS)= .399+.01tj*b 

GSIF(S) = .9975617977-. 018426966*S 

DIMENSION C3D(9no) »C3A(900) »KHO(900> »CDC0(9n0) »DACHD(900) »COFMD(90 
10) fSRD(900} »CDCA(900) .DACHA ( 900 ) fCDFMA ( 900 )* SRA (900 VrATIO (900 > 
DIMENSION VA0R(90(^) »VDOR(900) 

DIMENSION WVA(90U) »WVD(900) .CVA(90U) »CV0(900) *ALPHVA(900) »BETVA(90 
10) .BETVD(900) *ALPHVD(900) 

DIMENSION AT(900) .TEMP (900) »teT (900) rCT (900) .ALPHA! (900) .BETAK900) 
l.TOR(900) 

DIMENSION WD(9U0) 

DIMENSION ALT(yOO) »YA0R(900) »YDOR(900) 

DIMENSION WA(900) .CA( 900). ALPHAA( 900) »BETAA( 900) . 71(900 )» 12(900). 
173(900) .YA (900) » YD (900) .CD (900) .ALPHAD (900) .BETAD (900) .BD(gOO) .BA 
1(900) 

dimension AA(900) .RH0A(900) .CDA(900) .RHOD(900) .CDD(900) .AD(900) . 
17EMPA(900) .7EMPD(900) .VA(900) .VD(900) 

10 READ (5.302) (SA) 

302 FORMA! (13) 

READ(5.250) (C3.E) 

REAO(5.301) (MA.MD.N.AL7M) 

READ (5. 300) ( AA ( I ) .RHOA ( I ) . 7EMPA (I ) .CDA ( I ) . yA ( I ) . Ir 1 .mA ) 

READ(5.300) (AD(I) .RHOD(I) .lEMPO(I) .CDD(I) .VD(I) .1=1.MD) 
WR17E(6.251) (C3.E) 

WR17E(fa.304)MA»M0.N.AL7M 

WR17E (6.303) (AA (1 ) .RHOA( 1 ) .7EMPA(I) .CDA(1 ) . VA ( I ) . 1=1 .MA ) 

V»R17E (6.303) ( AD ( 1 ) »RHOU ( 1 ) »7EMPD( I ) . CDD ( I ) . VD ( I ) . 1 = 1 .MD) 

304 F0RMA7(///3HMA=. I4.3HMD=. I4.3H N=. I4.5HAL7M=.E10.5) 

303 FORMA! ( 4X .FO. 4 . 2X »E8 . 4 »2X . F8 ,4 .2X .F8 . 4. 2X .F9. 4/ ) 

GAMMA=1.4 

PI=3. 141592653 
C0NVSM=(GAMMA*.5)**.5 
DO flO 1=1. MA 
AT(I)=AA(I) 

80 !EMP(I)=!EMPa(I) 

DO 81 1=1. MD 

A!(MA4I)=AD(1) I 

81 !EMP(MA + 1)=!EMPDU) 

M!=MA+MD 

300 FORMA! (F7.3»Eb.3.F8.4.F5,3.F«.2) 

301 F0RMA!(3I3»F6.2) 

DO 70 1=1. MA 

70 WA(I)=1.0 

DO 71 1=1. MD 

71 WD(I)=1.0 
WR17E (6.305) SA 

305 FORMA! (4X. 17HS0DNDINb NUMBER =.I3/) 

V«R1TE(6.201) 

DO 2 1=1. MD 
YID=RHOD(I)*CDD(I) 

SXD=SXD+AD(I) 

V»Rl!t(6.100)AD(I) . YID 
Y=ALOG(RhOD(l)*CCD(I) ) 

YD(1)=Y 
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2 CONTINUE 

WRiTE(6»i:00) 

DO 1 I=1»NA 
Y1A=RH0A(I)*CDA(I) 

Y=AL0G(RH0A(1)*C0A(I) ) 

YAlDrY 

taRlTL(6»100)AA(I) >YIA 
1 CONTINUE 

WRITE <6»200) 

DO 60 I=I*MA 

60 WRITE (6*100) AA(i)»YA(I) 

WRITE (6*201) 

DO 61 I=l*MO 

61 WRITE (6*100) AD(I)*YD(1) 

N=bl 

DO 30 K=1*N 

30 ALT (K)=ALTM+.5*(K-1)410.0 
LT=0 
JT=0 
KFT=3 
KT=3 
KD=2 
KA=2 
KFA=2 
KFU=2 
KCA=2 
KCU=2 
LA=0 
LD=0 
JA=0 
JD=0 
KVA=3 
LVA=0 
LVU=0 
KVD=3 
JVA=0 
JVU=0 
KFVA=3 
KFVD=3 

CALL ORfHLS ( AT * TtMP* WT *MT*LT » JT*CT * ALPHAT*BETAT*KT * Tl *T2,T3» INDIT 

CALL ORThib ( AA * YA » WA »NA »LA *JA *CA * ALPHAA vBETAA »KA » Tl » T2* T3 * INDIA ) 
CALL ORTHLS (AD* YD*w 0*NU»L0*JD*C0*ALPHA0*BFTA0*KD*T1*T2*T3* INDIO) 
CALL ORTHLb ( AA » VA r M VA »MA »LVA » JVA *CVA * ALPHVA * BETVA *K VA r T1 * T2* T3 * IN 
IDIVA) 

CALL 0RTHLS(AD*VD»WVD*MD*LVn*jV0*CVD*ALPHVD*BETV0*KVD»Tl*T2*T3*IND 
IIVD) 

11=0 

WRITE (6* 103) IND1VA*1I*CVA( t) » (I1*CVA(1I<»1) *ALPHVA(1I)*BETVA(II)» 
11I=1*KVA) 

11=0 

WRITE (6*l03)lN0IVD*IlfCv0(l)*(II*CV0(IT-^l)*ALPHV0(Il>rBETVD(Il)* 
1II=1*KVD) 

11=0 

WRITE (6* 103) INDlA»iI*CA(l) » (II*CA(1 I-M) *ALPHAA(TI) *BETaA(IX) »TI=1* 
IKA) 

11 = 0 
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I 

MHlTk(6»103) lNDlD»iI»CO(l) f (XItCOai-^1) »ALPHAD(I1) *BETAD(II) »T1=1, 
IKD) 

11=0 

WRiTt(6»l03)INLllT»XI»CTa ) » ( 1 1 »CT ( II + l ) » ALPHA! ( II ) »BETAT ( 1 1 ) » 1 1=1 
IfKT) 

103 FORMAT(//I3///2X»l4r2XrF20.7f /(?X»I4»2X»3E2n,7) ) 

CALL COEFS <0A »CA » ALPHAA»BETAA»KCA»BA » T1 »T2» T3» IND2A ) 

CALL COEFS (OO »C0 » ALPHAD tBtTADfKCDjBD • T1 » T2 » T3 • IND2D ) 

CALL FITY(ALTiN» JA»CA»ALPHAA»BETAA»KFA»YA0R»T1»T2»IND3A) 

CALL F I T Y ( ALT » N » JL » CD » ALPHAD » BETAD » KFD » YDOR * T 1 » T2 » IN03D ) 

CALL FITY(ALT»N»JT»CT»ALPHAT»BETAT»KFT»T0R»T1»T?»IND3T) 

CALL FlTY(ALT*NrJVA»CVA>ALPHVA*PETVA>KFVArVA0R»TlfT2»lND3VA) 

CALL FITY(ALT»N»JVUrCVD»ALPHVD»BETVD»KFVD»VD0R»TlrT2»lND3VD) 

MRiTt (6r5U0) (K» ALT(K) *TOR(K) *VAOK(K) »VUOR(K) »K=1»N) 

500 F0RMAT(//3X»lHK»5X»6HALT<K>»8X»6HT0R(K)»aX»7HVA0R(K)aX»7HVD0H(K) »/ 
l(2Xtl3»Fl0.4*3E15.b)) 

V^RITE (6t450) ( I »UA ( 1 ) * 1 = 1 »KCA ) 

MRXTE (6»4bl) ( I »UD ( X ) » 1=1 »KCn) 

450 FORHAT <//4X» IhJ » IbX » 5HBA ( I ) »/ (2X r 15»3X»E20 .7) ) 

451 FORMAT ( //4X » IHl » IbX » 5HB0 ( X ) » / (?X » X5» 3X »E20 .7 ) ) 

WR1TE(6»202) 

DO 3 I=1»N 

EANM=24 .68+. 1235*AlT ( I ) -.000875*ALT ( I > *AlTI I ) 

R6=8314.34/EANM 

SRD( X)=VUOR( I)/(2.*R6*10R(I) )**,5 
IF (SPD(I)-1.25)lVb»195»X96 

195 S=bROU) 

CPU = S/tONVSM 
GSI = GSIF(CMD) 

CDFMD< X) =t2./(PX»*.5) )*(8./<3.#S)+fl.*S/15.-8.+S*S*S/2lO. )*GSI 
GO TO 197 

196 S=bRD(X) 

CMD = S/CONVSM 
GSX = GSXF(CmD) 

SI2=l./(b*S) 

Sl4=bI2*SI2 

CDFMD(X)=2.*(l.+Sl2-.2b*SI4)*GSI 

197 DALHD<X)=SRD(T)/C01\VSM 
DM=nACHD(I) 

COCA n)=.924.1b607l4/DM-. 3660714/ (DM*DM*DM) 

VRAT XO ( I ) =V AOR ( I ) / VDOR ( I ) 

SRA(i)=VAOH(l)/(2««RG«TOR(X) )««.5 
IF (SRA(X)“1.25) 9b»95»96 

95 S=bRA(X) 

CPA = S/tONVbP 
GSX = GSXF(CPA) 

CDFPA(l) =(2*/(PI**.5> )*(0./(3.»S)+e.*S/15.-8.*S*S*S/2lO.)*GSI 
GO TO 97 

96 S=SPA(I) 

CPA = S/CONVSP 
GSX = GSXFCCPA) 

SI2=l./(b*S) 

Sl4=SI2*bI2 

CDFMA(1)=2.*(1.+SX2-.25*SI4)*GSI 

97 OACHA(1)=SRA(I)/COnVSM 
DM=DACHAtI) 

CDCDd )=.92+,lb60/X4/DM-.36607l4/(DM*OM*DM) 
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RH0CUA=EXP(YA0R<1)) r 

RHOCDD=EXP ( YOOH Cl ) ) 

DENOMrCDL A m ♦CDPHU U ) -CDCD ( I > *CUFMA ( I ) 

RHO ( I ) = ( RHOCOA* ( CUR MU ( I ) -COCO <1 ) ) -RHOCDD* ( CDFMA ( I ) -COCA < I ) ) ) /OENOM 
RHOIN-l./RHOlI) 

XX=(RHOCDA*RH01N-COCA(1) )/(CDFPA(I)-CDCA(I) ) 

IF{XX)77«79»79 

77 C3A(I)=-1. 

GO TO 78 

79 X = XXF(CMA) 

C3A(I)=-(RHOIN**X )*AL06(XX) 

78 XY=(HH0CUD*RH01N-CUCD(1) )/(CDFMP(I)-CDCD(I) ) 

IF(XY) 87 » 89 f 89. 

87 C3D( 1 )=-l. 

60 TO 88 

89 X = XXF(CMU) 

C3U(I)=-(RH0IN**X )*ALOG(XY) 

88 AOVRD=RHOCUA/RHOCuU 
FMRA=CDFMA ( I ) /COFMlD < I ) 

CONRArCDCA ( I ) /CDCU i I ) 

TMilN=( (l.-A0VRU)/CA0VRU/(VU0R(I)**2,}-l./(VA0R(I)**2.) ) )/(2.«PG) 
MR1TE(6*102)ALT(1) »YAOR(l)*RHOCnA»YDOR(I) rRHOCOOrFMRAtAOVROrCONRA* 
ITMIN 

3 CONTINUE 

WRITE (6*333) ( ALT ( i ) * SRA < 1 ) r COCA ( 1 ) »COFMA ( I ) »OACHA ( 1 ) * 1 = 1 * N ) 

WRITE (6 » 334) ( ALT ( 1 ) » SRU ( I ) »COCP ( I ) »CDFMD ( I ) *DACHD ( I ) * 1=1 »N) 

WRITE (6* 33b) (ALTU) »C3A(1) »C30n> »VRATIO(I) rRHO(I) *I = 1*N) 

333 FORMAT (//4X»3HALT»9X»3HSRA*9X*4HCDCA»9Xf5HCDFMA»8X»5HDACHA*/(2X*Fl 
10.4*4E13.5) ) 

334 FORMAT ( //4X * 3HALT * lOX * 3HSRD » lOX » 4HCDC0 * 9X * 5HCDFMD » 8X » ShOACHD » / ( ?X » 
1F10.4»4E13.5) ) 

335 FORMAT (//4Xr3HALT»10X»3HC3A,lOX*3HC3Df8X*6HVRATIOr8Xr3HRHO»/(?Xf FI 
10.b*4E13.5) ) 

60 TO 10 

100 FORMAT(2X*F10.4t2X*E12.6) 

102 FORMAT (2X*F10.4*8El4.b) 

200 F0RMAT(///bX»5HAA(I) *10Xf 3HYIA//) 

201 F0KMAT(///5X»5HA0a) »lOX,3HYln//) 

202 FORMAT! ///fax »3H alt 1 12X » 2HYA » 8X *6HRH0CDA * lOX » 2HYP» lOX *6 hRH0C01) * 8X * 4 
1HFMRA»IOX»5HAOVRD*9X»5HCONRA) 

250 FORMAT(2E2U.10) 

251 FORMAT (//.4X»3HC3=»F20.10.4X*2HF=*E20, 10) 

END 


no 
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PROGRAM RHORAT 

This program takes falling sphere data and performs the usual 
analysis to determine the temperature and density profiles. The data 
In the region of overlap of the ascent and descent trajectories Is 
treated as If the atmosphere were Influenced by Isentroplc waves. 

Thus, differences In density at the same altitude will result In 
differences In temperature at that altitude according to the Isentroplc 
relations. 

The program also Inputs various standard atmospheres so that 
ratios of the measured density and temperature can be readily com- 
pared to the standard values. 


Ill 



OK UN, /TP Rh0RAT,UAHXXXXXXXXX,E:NVIR-DYN-GK,5i,l50 


XXF(S)» »39?^,016*S 
G5If(S» ■ .997S617977-.018H26966*S 
E Anhf ( 6 ) *29 •68+ • 1 235«G*«000S79«g*G 
parameter NPP-200 

DIMENSION MI5QO) , NAD (500) »NAA(BOO) ,TU(&00) »TU2(500) , THE (500) , 
1SIG(500) .PU(SOO) ,PU2(500) ,PME(500) ,siP<^UO) ,TMSR(500) ,SlSR(SOO) 
DIMENSION RHDD( 100) ,RAISO( 100) ,TRA( lOO) ,CISO( 100) 
dimension TTA(NPP) ,WTA(NPP) ,CTA(NPP) ,APTTA(NPP) ,BETTA(nPP) 
DIMENSION GR( 100) ,6RA( 100) ,ALPHG( I 00 > • BET AG ( 1 00 ) ,GRD( IQO) ,C6( 100) 
DIMENSION TINA(NPP) ,TIND(NPP) ,THA(NpP) ,THD(NPP) 

DIMENSION SRA(NPP) ,COFMA(NPP) ,0ACHA(NPP) , CDC A I NPP ) , SRO ( NPP ) .COFMD 
1 (NPP) ,DACHD(NPP) ,COCO(NPP) 

DIMENSION AA(NPP) ,RHOA(NPP) ,COA(NPP) ,RHOD(NPP) , COO (NPP) , AO (NPP) , 

1 TEMP A ( NPP ) , TEMPO ( NPP ) , VA ( NPP ) , VQ( NPp ) , 

1RH06A( 100) ,NAN( 100) ,RHORA( 100) , CdGA( 100) tCORA( 100) ,RH0GD( 100) ,N(}N 
1 ( 100 ) ,RHORD ( 100 ) ,CDGD( lOO) .COROMOO) 

DIMENSION HHLS(NPP) ,RHLSA(NPP) ,RHLSd(NPP) 
dimension RHOSA ( NPP ) ,RH0S0(NPP) 

DIMENSION ALS(NPP),RH0S(NPP) ,TES(NPp) ,Tl (NPP) ,T2(NPP) ,T3(NPP) ,FPS( 
INPP) .SMOL(NPP) ,W(NPP) ,CS(NPP) ,ALPS(nPP) iBETS(NPP) 

DIMENSION AT(NPP) 

1 ,XA(NPP) .XDINPP) ,C3A(NPP) ,C30(NPP) 

I ,RHA(NPP) ,CRA(NPP) ,APA(NPP) ,BEA(NPP) ,RH0(NPP) ,TEMA0( 100) 

1 ,TTA(NPP) ,TTO(NPP) 

CC ■ 307072000. •.69 
KTA ■ 8 
JTA ■ 0 
LTA ■ 0 

READ (5,601) (MS) 

READ (5,800) ( AtS ( 1 ) , RHOS ( 1 ) , TES ( 1 ) , FPS ( I ) , SMOL ( 1 ) , I - 1 , MS ) 

READ (5,600) (GR( I) , t>l ,MS) 

WRITE(6,599) 

599 FORMATllHO , 8X , ♦ ALT ♦ , 8X , * GR AV I TY * ) 

WRITE (6,601) ( ALS( 1 ) ,GR( I ) , !•! ,MS) 

600 FORMAT (10X,F10.5) 

WRITE (6,803) (MS) 

WRITE (6,802) ( ALS( I ) ,RHOS( I ) ,TCS( I ) ,FPS( I ) ,SM0L( I ) , I«I ,MS) 

801 format (110) 

800 FORMAT ( F 1 5 . 5 , E 1 5 • 5 , F 1 5 . 5 , E 1 5 . 5 , F 1 5 . 5 ) 

803 format ( ///,2X,3HMS«, I 10,2X, I3HU.S. STANDARD ,//, 9X , 3HAl5 , 1 2X , 9HRM0 
IS, 12X,3HTES, 12X,3HFPS, I 2X , 9HSM0L . // ) 

802 FORMAT ( 2 X , F 1 5 . 5 , E 1 5 . 5 , F 1 5 . 5 , E 1 5 • 5 , F 1 5 • 5 ) 

ER-.OOOl 

DO 60 I>) ,MS 

60 RHLS( I )«ALOG(RHOS( I ) ) 

KKK«6 

CALL ORTHLS (ALS,RHLS,W,MS,0,0,CS,ALPS,BETS,K)CK,T1 *T2,T3,IN01 ) 
GAMMA* 1 . 9 

C0NVSM>(GAMMA«.5)*«.5 
PI«3. 191592653 
DO 89 Ia80,120 
M( I )»0 
TU( I )*0. 

TU2(I)aO.O 
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T M t ( I » « 0 . 0 
5 1 Cl ( I ) sU« 
pj ( n *0. 
pu<! ( I ) »n, 

PME ( I ) «0. 

0M 5 I P ( I ) -0 I 

10 KLAO (B,302) ISA) 

ROAD (5.310) ( HAD ) 

ROAD (b,2UQ) IC3.C) 

R0A0(b,301 ) IMA.MD.N.ALTM) 

R0AD(S.300)(AA(I),RH0A(n»TEMPAII),cDA(I),VAIl),I«l,MA) 
ROAO|S,300)(AOI1),RHOD|I) , TEMPO! 1) , CDD ( I ) , VD ( I ) , 1 » I , MO ) 
tea ■ TEMPA I II 
TOD ■ TEMPO ( 1 ) 

00 M8 I ■ 1 ,MA 

NAA ( 1 ) >AA ( I ) ■^•2 
4B AT ( I ) > AA < I ) 

DU 85 1 « 1 |MD 

TOMAD ( I ) -0.0 
RHDO ( 1 ) >0*0 
as NAD! 1 )»A0( I )4^«2 
DU H9 1 a 1,5 
•49 WTA I I ) ■ .5 

DO <47 I ■ 6.MA 
•47 »^rA( I ) a I .0 

CALL ORTHLS I AL S , 6R , W , MS , 0 , 0 , C G , A LP hG , BE T A G i 4 , T I , T 2 . T 3 , ! ND 1 ) 
CALL FITY ( AA ,MA ,0 ,CG , ALPHG.BETAG ,4 ,6RA ,T I , T2 , IN03 ) 

CALL FITY ( AD,MD,0,CG,ALPHG,BETAG,4,GRD,TI ,T2, 1N03) 

6U 1 FOHMA T I 4X , 2F I 0.5 ) 

WRITE 16.601) ( AA I 1 ) .GRA! I ) • 1-1 .MA ) 

WRITE (6,601) I AD I I ) ,GKD( I ) . lal |M0 ) 

DU 5 I a I , 1 00 
CORA ( I ) aO. 0 
5 CURD I I I »0.0 

WRITE (6,305) SA 
WRITE(6,311)|MAD) 

WRITE 16,201) IC3.E) 

WR1TF(6,304)MA,M0,N,ALTM 

WR|TF(6,303)(AA(I),RH0A(I),TEMPA(I).CDA(I),VA(I),Ia|,HA) 
wRITE(6, 303)140(1), RH0D(I),TEMPD(I),C 00 (l),VD(l),Ial, MO) 

305 FORMAT I4X, 17HS0UNO1NG NUMBER a,I3/» 

30 1 FORMAT ( 3 I 3 ,F6,2 ) 

300 FORMAT I F 7 . 3 , 0 8 . 4 , F 8 . 4 , F 5 • 3 , F 8 . 2 ) 

311 FORMAT ( 4X ,4HMAD- , 1 4/ ) 

310 FORMAT (12) 

304 format ( ///3HMA« , I 4 , 3HM0« , 1 4 , 3H N ■ , I 4 , &H A L T M» , E 1 0 • & ) 

303 FORMAT ( 4 X , F 6 . 2 , 2 X , E 7 . 2 , 2 X , F H , 0 , 2 X , F 5 . 3 , 2 X , F 6 . I / ) 

302 FORMAT (13) 

200 FORMAT (2E20.10) 

201 format ( // , 4X , 3HC3- ,E2U. 1 0 , 4x , 2HEa ,£2 o. 1 0 ) 

CALL FITY(AA,MA,0,CS,ALPS,BETS,KKK,RHL5A»T1,T2,1ND3) 
call FITY ( AD ,M0 ,0 ,CS , ALPS ,BET5 ,KKK .RHLSO , T 1 ,T2 , I Nd4 ) 

MNO-O 

DO 55 I a 1 ,MD 
55 RHOGD ( I ) aRHOO ( I ) 

76 DO 12 I a 1 ,MA 
ALT AaA A ( t ) 
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R6-831 **-. 3M/EANHf ( aLTA ) 

SRA( I )»VA( I i/( 2 ««rg*tehpa( I n**«s 

IF (SRA(l)-}.2&) 95i95|96 
9b S»SRA( I ) 

CMa ■ S/CONVSM 
ftSi ■ gsificma ) 

CUFMAi 1 ) ■(2,/(PU*.5n*(a,/{3,*S)4.p,*S/t&,.e»«S*S*S/210, )«6S1 
GO TO 97 

96 S»SRAm 

CMA « S/CONVSH 
QSI ■ GSIF(CMA) 

SI2«i ,/(S«S) 

SIH-S12»SI2 

CUFmA( I )>2t*( t •♦Sl2-.2b«SI*< )*GSI 

97 dacha J 1 ) «SRA ( 1 ) /CONVSM 
DM«OACHA ( I ) 

CDCA(I)b* 92«« 1660 7 14/0t1-«36607 1M/|DM*0M*DM) 

DC«C0FMA( I )«CDCA( I I 
NN«20 

RHOO«RHOA { I I 
K«-I 

XA I 1 ) ■ XXF (CMA ) 

X ■ XXFICMAJ 

C3AII) ■ ,212*(CC/E;aNMF( ALTAM**XA( I ) 

C3 ■ C3A ( I ) 

30 FX»RHOA« I »*COA( n/KDCA( I »*OC*EXPI-c3*«RHOO»*X )») 

CALL WEGIT(RHOU,FX|E(K,NN) 

GO TO (30,31 ,32,33) ,K 

31 RHOGA( I )-RHOO 
GO TO 3H 

32 RHOGA ( I ) *0*0 
GO TO 3H 

33 RHOGA ( I ) «- 1 .0 
39 NAN ( 1 )'«NN 

RH0RA(I)«RH0GA(I)/RH0A(I) 

CD6A ( 1 ) -COCA U )*DC*EXP(-C3*(RH0GA( I )»*X ) ) 

RHOSA( 1 )>EXP(RHLSA( I ) ) 

CORA ( I ) -RHOGA ( 1 ) /RHOSA ( I ) 

12 continue 

IFlMAD-BO) 20,21,21 

20 DO S8 JK-1 ,MA0 

TEMAO ( JK ) -TEMPA ( MA»MAO-f JK ) 

RHDD( JK ) -RHOGA (MA-MA0<»JK ) 

5B TEMPD(JK)-TEMAD(JK) • ( RHOGD ( JK ) /RHOO ( JK ))•*(•*« ) 

GO TO 22 

21 TEMPD( 1 ) - TED 

22 DO II I - 1 ,HD 
ALTO-AO ( 1 ) 

RG-B3|S.3‘(/CANHF(ALTD) 

SRDI I )-VU( 1 )/I2.»RG»TEmP 0( I 1 )**,5 
IF |SRD( 1 )-l .25) 195, 19b, 196 

1 95 5»SR0 ( I ) 

CMD - S/CONVSM 
GSl - GSIF(CMD) 

CDFMD(l) -(2,/lPI««.b) )*(8./(3.»S)«B**S/lS."e.*S*S*S/2)0« )*651 
GO TO 197 

196 S-5R0( 1 ) 
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CMD ■ S/cONVSM 
6bl ■ GSIF(CMUJ 
SI2-1 ,/(S*S) 

SIH»SI2*SI2 

CDFMO( I )»2.*M .♦SI2-.2b*SIH)*6SI 
197 DAChD( I )«SRO ( I ) /CONVSM 
OM-OACHD « j r ■■ ■ 

COCD( I |«.92«« 16607 1H/DM.,366071‘i/(DH«DM«DH) 

DC»COFMO( I l-COCO( I) 

NN-20 

K»-l . : 

RH0G«KH0D( I ) 

XD ( I ) ■ XXF ( CMD ) 

X ■ XXF (CMD » 

C3D(1) ■ .212*(CC/EANMF( ALTD) »»*XD( I ) 

C3 ■ C3D( I) 

MO FX*RH00 ( I ) "CDD ( I) / ( COCO ( I (♦DC^EXP ( -C3» < )M 

CAUL WEG I T ( RHOU tFX •£ |K tNN ) 

GO TO ( 40. HI .H2.H3) ,K 
H 1 RHOGD ( I ) «RHOO 
60 TO 44 

42 RHOGD ( 1 ) «0.0 
GO TO 44 

43 RHOGO( n«-l *0 

44 NON ( I ) «NN 

RHORO ( I ) ■RHOGD ( I ) /RHOD ( I) 

C060 ( I ) -COCD ( I) ■i.DC*EXP ( -C3* { RH06D( I ) ••X )) 

RH0SD( I )«EXP(RHLSD( I ) ) 

CORO { I ) -RHOGD ( I ) /RHQSD ( I ) 

RAIS0(I)a(RH0D(l ) /RHOGD ( 1) )»*,H 
II CONTINUE 

MRITE (6,902) 

RRITE (6,900) ( AA I I ) ,RhOGA ( I ) .RhORA ( I ) , COGA ( I ) , COR A ( { ) , N AN ( I ) , 

I XA( I ) ,C3A( I I . )■) ,MA) 

WRITE (6,901 ) 

WRITE (6,900) ( A0( I ) ,RhOGO( I ) ,RhDD( I ) , CDgD ( I ) , COrO ( I ) , NON ( I I , 

1 X0( I ) ,C3D( I) , I«l ,M0) 

BB-0.0 

DO 6) X«2,MA 

BU-BR.(AA(K)-AA(K.ii)*(RHOGA(K)»GRA(K)>RHOGA(K. |)»GRA(K.t)| /AL06 
1 (RH06A(K )«GRA(K )/(RH0GA(K-1 )*GRa(K-| ) ) ) 

61 T1NA(K )»BB 
BB-0,0 

DO 63 K«2,M0 ' 

BB«B8-(AD(K)-A0(K-1 ) )»(RH060(K)*GRD(K)-RH0GD(K-1 )*GRD(k- 1 ) I/AUOG 
2(RH0GD(K )*GRD(K )/(RHOGD(K-t )«6RD(K-| ) ) ) 

63 T I NO ( K ) -BB 

THA ( I ) -TEMPA ( 1 ) 

DO 64 KK«2,MA 
ALTA«AA(I^K) 

RG«8 3 I 4, 34/EANMF ( ALTA ) 

64 THA(kk)»TINA(KK)»100 0./(RH0GA(kk)*RGI ♦RHoGA ( 1 ) *THA ( 1 ) /rHOGA ( kk ) 
THD ( I ) -TEMPO ( ) ) 

DO 65 KK-2,HD 

- ALTO-AD ( XK )‘ ■ - 

RG-8314.34/EANMFI ALTO) 

65 THD(KK)«TIND(XK)*lOaO./(RH06D(KK)*RG)*RMOGO( 1 )*TH0( I )/RHOQD(KK) 


115 


DO 57 1 « J ,mO 

TRAC I )«TE;MAOn »/TEhPD«l ) 

57 CISO( 1 )«RAISO( I )/TRA( I ) 

WRITE (AfAlO) < AA< I > »TEHPa< I > »Tha( I) » !■» iHa> 

WRiTE(6,6l I M AOC I » ,TE mPDC I » iThOC i ) iRAiSOC i ) .TEHAOC j } ,C|SOl I 

I ) , I»l ,MO). 

DO 7 t I > I ,MA 

0|f*«ABS( (TEMPAC I )-Th*< I » »/ThM I ) )*DIF* 

71 TEMPAC I )»ABSCTHACI » ) 

IF (HAD-5Q) 23,2H,24 

23 MA1*MaO-»1 
GO TO 25 

24 MAI > 1 

2G DO 72 I«MAi,MO 

DIFD-ABSC C TEMPO! I )-THDC U l/THDC 1 M^OlPO 

72 TEMPDC I »«ABS(THDC I> ) 
sumo»oifa*oifo 

TEMPA C I ) ■ TEA 
IF (SUMD-ER) TR.TR.TS 
75 MNO-MNO+I 

IF (MNO-5) 78,78,74 

78 GO TO 76 
7H CONTIMUE 

DO 8 1 I > 1 ,MD 
JC>NAD ( I ) 

M( JC)«MC JC|4l 

TU2C JO »TU2CJC )♦ TEMPO ( tf»2 
PU2(JC)»PU2(JC)>CDR0(I)««2 
TUC JC)«TU( JO+TEMPOC I ) 

PUC JO-PUC JO*CDRD( n 

81 CONTINUE 

DO 62 I >1 ,MA 
JC-NAACn 
M( JO-MCJC)*! 

PU2CJC)»PU2CJ0^CDRACI)«»2 
TUC JC)»TUCJC)^TEMPA( I ) 

PU ( JC I -PU ( JC ) *CORA ( I » 

TU2C jC)«TU2(JO*TEMP AC n«*2 

82 CONTINUE 

00 83 I«80,120 
TMEC I )«TUC I )/MC I ) 

SIGCl)»CCTU2Cll/MCl))-TMEm**2J««.5 
PME c I )«PU( I )/M( I » 

J ■ 1-79 

TMSRC I ) » TMEC 1 I/TESC Jl 
SISRC I ) * SIGC i )/TESC JC 
8 3 5IPC I »«C CPUZC I »/M( n J-PME ( n*»2»*«.5 

WRITE! 4, 650 I C I ,Ml 1 ) ,TME< I> ,SIG< I > ,PmE< I > iSlPC I ) ,TMSRI I ) ,S1SR( I ) , 

II >80,120) 

650 format I 6X , IHI ,6X ,4HMC I I ,7X,6hTME C I) . 9X,6hSIG(P» 9X,6hPME(D, 

1 9x ,6HS1PC I ) ,9X ,7 hTmSRC I ) ,9X,7HSISRc I ) / / ( 5X , I 3 , 5X , I 3 , 6E 1 5 . 7 )) 

6 10 format C 6X ,5HAA C I ) , I OX , 8HTEMPA ( ] ) , mx,6HTHA( 1 I // ( 4X i F 1 0 • 5 , 2E20 1 1 0 ) 
3 ) 

61 I F0RMaT(6X,5HA0< I ) , 1 OX , 8HTEMPD ( I ) , 1 3 X » SHTEMAD I I 1 t9X,8HRAlS0l I > t 
t 12X ,6HTRAC I ) , lUX ,7HCIS0C n//(4X,F10.5.5E20t 10) ) 

900 format I 2X ,F 1 0*H , 4E 1 9.6 , I I0,2E 1 <I*6 1 

901 FORMAT C /// ,6X , 3HALT ,9X ,SHRH0GDt9X ,5HRHDD tlOX, 4HC0GD , 1 OX , 4HCDRD t 
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1 1 1X,3HADN«9X**XA* t 13X 

,«C3a* 

1 




902 

exQT 

FORMAT (///|8X(3HALT« 
1 1 IX ,3HAAN ,9X . *X0* , 13X 
aO TO 10 
END 

91 

9X*5HRH0<>A.9X( 
, 'CSD* 1 

5HRH0RA. lOX 

.9HCDGA. 10X.9HC0RA 

eo* 

1.999 

-E5180 

65 

9.065 

-E328 

969 

81 . 

1.662 

-E5180 

65 

9.888 

-E328 

969 

82. 

1.382 

-E5180 

65 

5.877 

-E328 

969 

83. 

1 .150 

-E5180 

65 

7*067 

-E328 

969 

8M. 

9.563 

-E6180 

65 

8.996 

-E328 

969 

8S. 

7.955 

-E6180 

65 

1 *021 

-E228 

969 

8&. 

6.617 

-C61B0 

65 

1.228 

-E 228 

969 

87. 

5. 509 

-E6180 

65 

1.976 

-F 228 

969 

88. 

9.579 

-E6180 

65 

1.779 

-E228 

969 

89. 

3.819 

-E6180 

65 

2. 133 

-E228 

969 

90. 

3.170 

-E6180 

65 

2.563 

-E228 

96 

91. 

2.598 

-E6183 

63 

3.127 

-E228 

96 

92. 

2.137 

-E6186 

62 

3.802 

-E228 

96 

93. 

1*763 

-E6189 

59 

9.607 

-E228 

96 

9M. 

1.959 

-E6192 

56 

5.566 

-E828 

95 

95. 

1.211 

-E6195 

51 

6.702 

-F228 

99 

96. 

1 *008 

-E6198 

95 

8.052 

-E828 

99 

97. 

8.915 

-E7201 

37 

9.693 

-E828 

92 

98. 

7.099 

-E7209 

28 

1.151 

-E128 

91 

99. 

5.911 

-E7207 

1 6 

1.371 

-F128 

90 

100. 

9.979 

-E7210 

02 

1.629 

'E128 

88 

101 • 

9.169 

-E7219 

86 

1.996 

-E128 

86 

102. 

3.993 

-E7219 

66 

2.316 

-E128 

83 

103* 

2.995 

-E7229 

93 

2.799 

-E128 

81 

109. 

2.992 

-E7229 

18 

3.290 

-E128 

78 

105. 

2.117 

-E7233 

90 

3.810 

-E128 

75 

108* 

1.809 

-E7238 

58 

9.965 

-E128 

72 

107* 

1*593 

-E7293 

23 

5.215 

-E128 

68 

108* 

1*323 

-E7297 

85 

6.07 1 

-E128 

669 

109* 

1.139 

-E72S2 

99 

7.095 

-E128 

60 

110* 

9*829 

-E8257 

00 

8.150 

-E128 

56 

111* 

8*360 

•E8266 

99 

9.568 

-E128 

51 

112* 

7.153 

-E8276 

85 

1.117 

28 

97 

113* 

6. 153 

-E8285 

20 

1.296 

28 

92 

119. 

5.321 

-E8299 

92 

1.996 

28 

37 

115* 

9.623 

-E8303 

78 

1.719 

28 

32 

116* 

9.035 

-E8313 

01 

1.966 

28 

27 

1 17. 

3.536 

-E8322 

19 

2.239 

28 

22 

118* 

3.112 

-E8331 

33 

2.590 

28 

17 

1 19. 

2.798 

•E8390 

93 

2*870 

28 

12 

120* 
80* 
81 . 
82. 
83. 
8H* 
85. 
86 * 

87. 

88. 

2.936 

9,569 

9.561 

9.558 

9.555 

9.552 

9.550 

9.597 

9.599 

9.591 

-E8399 

99 

3.233 

28 

07 
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0UU3M 1 

QUO 

0003H2 

000 

0003H3 

000 

0U03M4 

GOO 

0U034b 

000 

0UU346 

000 

000347 

000 

000348 

000 

000349 

000 

000350 

000 

00035 1 

000 

000352 

000 

000353 

000 

000354 

000 

000355 

000 

000356 

000 

000357 

000 

000358 

000 

000359 

000 

000360 

000 

00036 1 

000 

000362 

000 

000363 

000 

000364 

000 

000365 

000 

000366 

000 

000367 

000 

000368 

000 

000369 

000 

000370 

000 

000371 

000 

000372 

w u u 


89. 

9.538 

90* 

9.535 

91 * 

9,532 

92. 

9,529 

93. 

9,526 

94. 

9,523 

95. 

9.520 

96. 

9,517 

97. 

9,514 

98. 

9,511 

99. 

9.508 

100. 

9,505 

10 1. 

9.502 

102. 

9,499 

1 03. 

9,496 

1 04. 

9,493 

105. 

9,490 

106. 

9,488 

107. 

9,485 

106. 

9,482 

1 09. 

9,479 

1 10. 

9,476 

111. 

9.473 

1 12. 

9,470 

1 13. 

9,467 

114. 

9,464 

1 15. 

9.461 

1 1 6 . 

9,458 

117. 

9,455 

1 18. 

9,462 

119. 

9,449 

1 yn . 

9.447 


118 



PROGRAM AFILIP-HIGH CM 


Thle program computes values of drag coefficient for various 
values of Knudsen number and Reynolds numbers for Mach numbers above 1 . 
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yHUN,/TPC AFILlPtUAHXXXXXXXXX,OR8IT,5,HOO 


HIGH CM 

SFOR.IS MAIN. main 

6SIF(S) « . 9V7&6I7977-.0ia‘»26966«S 
XF(S) ■ .399+.016*S 

CDFMLF(5)*I2»/<PI«**5» )*<8./(3.*SI+fl.*S/I5.-8.*S»S*S/210.I*GSI 
COFMHF(S) • 2.« ( I . + l ./ I S»S >-.25/ « S*S*5*S ) )*GSI 
COCF(S)». 92*. 1660719 /S-. 3660719 /(S*S«S) 

DIMENSION XH( 100) ,RE( 100) «C0( 100) ,BKN( 100) I 
1 COLS( 100) ,CR( 100) 

J»1 » 3.19159265 
GAMMA * 1.9 

CONVSM » (GAMMA*. 5)»*.b 

10 READ(5. 1 .END«100)CM,NI 

1 F0RMAT(F10.0» 15) 

R£AD(5,33) Xt.EL.GSI 

33 FORMAT ( 3F 15.0 ) 

RLA0(5.2) (XM(I),RE(1).C0(I).I«1.NI) 

2 format ( 3F20.0 ) 

GSI « GSIF(CM ) 

SELK « n.o 
SELK2 » 0.0 

SY ■ n.o 
SYKN > Q.O 
WHl TE ( 6 , 1 1 ) 

11 FORMAT! iHl . 13X , 'PHESENT EXPERIMENTAL SPHERE DATA*) 

WRITE(6,5)CM 

5 FORMAT! lH0,/,/|25X. *CM > • * F6 , 9 , / , / , 1 I X , • XM • . I 7 X , * RE » , 1 9X , ♦ CD » » / ) 
WRI TE ( 6 .6 ) ( XM( 1 ) ,RE ( 1 ) ,C0( I ) » I«1 ,NI ) 

6 F0RMAT(5XiF10.9,10X,F10.9,I0X,F10.9) 

AGbGAHMA** ( -.5 ) 

7 DO 20 I «1 ,NI 

A » .999*18. /PI )**.5 
A ■ I ./A 

BKN( I ) > XM( I )/(RE( 1 )*AG) 

BKN ( I ) « A*BKN ( I ) 

SR ■ XM ( I ) *CONVSM 
IFlSR-1.25) 95,95,96 

95 CDFM * CDFMLF(SR) 

GO TO 97 

96 CDFM B CDFMHF(SR) 

97 6M B XM ( I ) 

CDCN B CDCF(GM) 

DC « CDFM-COCN 
UC > C0( I )-CDCN 
CK( I ) » UC/DC 
WRITE(6,9)0C,UC,CR( I ) 

9 format (/,/, 1 H ,«DC •» ,E20. 10,5X, *UC b . , E 20 . 1 0 , 5 X , ♦ C R ( I ) >*,E20.I0) 
20 CONTINUE 
E * .212 
DO 91 JB] , 101 
X«.005*(J-1) *0.35 
FXbDSUX(CR,BKN,E,N1 ,X) 
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I 


F X-FX-X 

SJ WK 1 TE ( 6 , 1 2 ) X ,FX ,NI 
X « XF(CM) 

00 H2 J« 1 t 1 01 

E».00S*< J-1 )♦. I 
F£«DSUE ( CH iBKN ,E iNl ,X ) 

FE-FE-E 

42 WR1TE(&, 13)E,FE,NI 
X » XF ( CM ) 

E - . 2 1 2 

12 FORMAT!/, IH ,*X ■ ♦ , E 20 . I 0 , 5 X F X « •', E 20 • I 0 » 5 X , * N ■ »,13» 

13 FORMAT!/, Iri ,»E ■ • , E 20 . 1 0 , 5 X , > F E - t , E 20 . 1 0 , 5 X , • N . *,13) 

SW0»0,0 

DO 7S I « 1 ,N1 
5R«CM*C0NVSM 
GM«CM 

1 F ! SR- 1 . 2b 1 SB • ,86 

8b CDLSll) ■ CDCF ! GM ) ♦ ( ( COFMLF ! SR ) -COCF ! GM > > ‘EXP I -E/( BKN ( I ) ) t-X ) ) 

GO TO 70 

86 COLSII) ■ CDCF !GM)*(C0FMHF1SR.)-C0CF|6M» )*EXP1-E/(BK:N( I ) »»»X) 

70 S^JD»SOD^(CDLS! 1 )-CD( I) )**2. 

7b CONTINUE 

RMS» ( SQO/N I I •• .5 
WR I TE I 6 , 7*n 

74 F0RMAT!/,/,/,9X,»C0LS»,18X,*C0',17X.*BKN*,/) 

WH1TE(6,73) (CDLSINl ,CD!N) ,8KN!N) ,N>1 ,NI ) 

73 FORMAT ! l.H ,3E2Q, 10 1 

WR 1 TF ( 6 , 3 ) CM ,N 1 , X ,E , 6S I 

3 FORMAT!/,/,/,/./, IH ,9X,*CM»,l8X,*Nl*,/,/,lH ,7X,F6.4,16X,I2,/,/,/ 
1,1H ,9X , »X» , IVX , *E» , 18X , »GSI «,/,/, IH ,3£20,10) 

WK1TE!6,71 IRMS 

71 FORMAT !/,/,/, ♦ RMS* ♦ , E20 . 1 0 ) 

S I « CM 

S2 ■ CM*C0NVSM 
DCF ■ CDCF ! S 1 I 
DFML - CDFMLF I 52 1 
OFMH ■ CDFMHF1S2) 

WR1TE|6,76) OCF , OFML, DFMH, 51, S2 
76 FORMAT !/,/,/, IM ,♦ CDCF* ♦ ,E20. 1 0 » / * / * 1 H ,♦ C DFMLF* ♦ , E 20 , 1 0 , 

I /,/,lH ,» CDFMHF** ,E20. 10,/,/, IH ,* S 1 ■ ♦ , E 20 , 1 0 , / , / , 1 H , 

1 - * S2»* ,E20. 10) 

DO b3 L* 1 ,B 
00 BO K*l,9 

BKNIK) ■ !0.01-»0*0l*!K-l ) )•! 10**«(L«i ) 1 
BKNIK) * BKN ! K ) • ( 1 0, •• ! -3 1 ) 

SR « CM*CONVSM 
GM ■ CM 

IFISR-1 .25)55,55,56 

55 COLSIK) » CDCFIGM)*! (C0FMLF!SR)-C0CF(GM) )*EXP(-E/I BKN!K))**X)) 

GO TO 51 , 

56 CDLSdC) • CDCF1GM|+(CDFMHF!SR» -CDCF(6M) )*EXP(-E/1 BKN(K))*»X) 

51 WK1TE!6,52)CM,SR,BKN(K) , COLSIK) 

52 FORMAT ( 1 HO , 1 6X , »CM* , 1 4X , *SR ♦,/,/, IH , 1 2 X , F 8 . 4 , 8 X , F 1 0 , 5 , / , / , / , 1 H , 

1 9X , •BKN* , 1 7X , «COLS« ,/ ,/ , IH , 3 X , E 1 5 , 8 , 5 X , E I 5 . 8 ) 

50 CONTINUE 

53 CONTINUE 
GO TO 10 
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100 STOP 

rUNCTlON OSUX(CRfBKN,E|NI ,X ) 
dimension CR( IOO) |BKN( 100) 

SX ■ 0.0 
DO 65 I ■ I |NI 
BKNX ■ BKN( I )*«X 

SX ■ SX*(CR( I l-EXP(-E/8 XNX) )*CEXP|-e/8KNX) )*(aL 06C 
6b CONTINOE 

OSOX ■ X>SX 
RETURN 

rUNCriON OSUEICRiBKN,E«NI ,X> 

DIMENSION CR( 100) ,BKN( 100) 

SE ■ OtO 
DO 66 I >1 ,N1 
BKNX ■ BKN( I )**X 

SE • SE^ICRI I )-EXP(-E/BKNX) )*(EXPI-E/BKNX) )/BKNX 
66 CONTINUE 

DSUE ■ E^SE 

RETURN 

END 


8KN( I ) ) )/BKNX 
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PROGRAM AFILIP-LOW CM 


This program computes those same values as the AFILIP-HIGH CM 
except that only Mach numbers lower than unity are employed. 



9KUN,/TPC AF ILIP.UAHXXXXXXXXX ,0RBIT,H,300 


LOW CM 

UFORiIS MAlMiMAlN 

XF(S) ■ .399>t0l6«S 

GSIFISJ ■ .9975617977-. 01892A96A*S 

C0FMLF(S)»(2./(PI**.5)»*<8./13.*S»-S.»S/l5.-e,*S*S*S/2l0.)*6S| 
COFMHF(S) ■ 2.*l 1 .♦! ./lS«S»-.25/(S*S*S*SM*GSl 
COCF(S) ■ 0.M0297-0. 1H2H799*S*S-0.M5950669*(S**R. ) 

DIMENSION XM(100),RE(100).CD(I00).BkN(100)i 
1 COLS(IOQ) «CR( loot 

PI ■ 3.19159265 
GAMMA ■ I.M 

CONVSM » (GAMMA*. 5)**. 5 

10 REA0(5i I .END«I00»CM,NI 

1 FORMAT ( F 1 0.0 • 15 I 
RCA0(5.33) XL.ELiGSl 

33 F0RMAT(3F15.0» 

REA0(5.2) (XM(I).RE(I)tCD(I),I«I|Nn 

2 FORMAT ( 3F20.0 J 
6S1 » 6SIF(CMJ 
SELK « 0.0 
SELK2 ■ 0.0 

SY - J.O 
SYKN « 0.0 
MR I TE ( 6 , I I ) 

11 FORMAT! IHI ,13X, ‘PRESENT EXPERIMENTAL SPHERE DATAM 
WRl TE(6,5)CM 

b FORMAT ( I HO , / I / , 25X I ‘CM ■ • » F 6 . 9 , / , / , I I X . • X M » , 1 7 X , • R E » , 1 9 X , • CO • . / ) 
WRITE(6,6)(XM(n,RE(l>,CD(I».I*l,NI) 

6 FORMAT(5X,F10.9,ipX,F10.H»IOX,Fl0.9) 

AG-G AMMA»« ( 0.5 ) 

7 DO 20 I>I ,NI 

A » .M99*(8./PI»**.5 
A ■ I ,/A 

BKN( I I • XM( I »/(RE( I )*AG1 
BKN ( I ) « A*BKN ( I ) 

SR ■ XM ( I ) *C0NVSM 
IF(SR-1.2S) 95,95,96 

95 CDFM » CDFMLF(SR) 

GO TO 97 

96 CDFM « CDFMHF(SR) 

97 GM » XM ( I ) 

CDCN » CDCF(GM) 

DC ■ CDFM-CDCN 
OC ■ CD ( I ) -CDCN 
CR ( I) » UC/OC 
WR1TE(6,9)DC,UC,CR( I ) 

H format «/,/. IH ,‘OC -* ,E20. I0,5X , ‘UC ■ • , E 20 . I 0 , 5 X , ♦ C R ( I ) ■*,E20.10 
20 CONTINUE 
E ■ .212 
00 9 1 J« 1 , I 01 

X-.005* ( J-1 ) *0.35 
FX*0SUX(CR«BKN,E,NI,X) 
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FX»FX-X 

HI MKITE(6,12)X,FX,NI 
X « XF(CMJ 
DO H2 J«1 , 101 
E«.005*< J-1 1 

FE«0SUE(CH,BKN,E,NI ,X» 

Ft»FE-E 

H2 WRITE(6. 13)E»FE|N1 
X ■ XF<CM) 

E > f2l2 

12 format (/, IH ,* X ■ * tE20. 1 0 ,5X , *FX ■ » , E 20 . 1 0 , 5 X , • N « *|13) 

13 F0RMAH/|1H ,»E ■ • i E20 . 1 0 , 5 X » ♦ FE ■ ♦ , E20 . 1 0 » 5X , * N ■ *,I3) 
soD«n,o 

00 75 

SR«CM*CONVSM 

(jM»CM 

IF(SR-1 .25)85i65»86 

8b COUS(I) ■ CDCF(GM)^( 1CDFMLF(SR»-C0CF«GMM*EXP(-E/IBKN( I ) )»«X) ) 

GO TO 70 

86 COLS(l) ■ CDCF(GM)*(COFMHF(SR J-COCF«GM) )*EXP(-E/(BKN( I ) )«*X) 

70 SQD«SQD*(COUS( M-CO( I n**2. 

7b continue 

RMS»(SQD/Nl )«*.b 
WR I TE ( 6 ,7H ) 

7H F0RMAT(/»/»/t9X,*C0LS*il8X,*CD*il7X.»BKN*,/) 
WRITE(6,73)(C0LS(N)|CD(N)|BKN|N)tN>l,Nl) 

73 FORMATIIH ,3E20.10> 

mH 1 TE ( 6 , 3 KB |NI ,E ,GSI 

3 FORMAT(/i/./i/./i IH ,9X, »CM» , 18X, »NI »././• IH , 7 X , F 6 . H , 1 6 X t 1 2 i / , / , / 

1,1H ,9X , *X» , 19X , *E» , 18X , 'GSl *,/,/. IH ,3E20tl0) 

WKITE(A,71 IRMS 

71 FORMAT I /»/ t /• * RMS»» |E20. 10) 

5 1 ■ CM 

52 ■ CM*CONVSM 
DCF ■ CDCF ISl ) 

DFML ■ CDFMLF(S2) 

OFMH ■ C0FMHF(S2) 

WRITE(6,76) DCF, OFML, OFMH, SI, S2 
76 FORMAT (/,/,/» IM ,♦ C DC F - ' , E 20 . 1 0 » / • t 1 H ,• COFMLF- ♦ , E 20 . I 0 J 

1 /i/,lH COFMHF** ,E20. I 0 , IH ,♦ S I ■ • » E 20 , I 0 , / , / » 1 H , 

1 • 52«* ,E20. 10) 

00 53 L-1 ,5 
00 SO K«1 .9 

BKN(K) « ( 0*0 1 -^0 *0 1 * ( K- 1 ) ) * I 1 0*«« (L- 1) ) 

BKN(K) ■ BKN(KI«( 10,««(-3M 

5R ■ CM*CONVSM 
GM ■ CM 

IFISR-l .25)55,55,56 

55 COLS(K) ■ COCFlGM)+( (COFMLF(SR)-COCF(GM) )*EXP(-E/( BKN<K))«*X)) 

GO TO 51 

56 COLS(K) « CDCF(GM)+ICDFMHF(SR)-CDCFIGM) )*EXPI-E/( BXN(K))**X) 

51 WHITE(8,52)CM,SR,BXN(K),C0LS«K) 

52 FORMAT I IHO , 16X , 'CM* , 1 HX , 'SR* ,/ ,/ , IH , I 2 X , F 8 . H , 8 X , F 1 0 . 5 , / , / , / t 1 H , 

1 9X,*flXN»,17X,»C0L5»,/»/tlH , 3X * E 1 5 , 8 , 5X , £ 1 5 . 8 ) 

50 CONTINUE 

53 CONTINUE 
GO TO 10 
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100 STOP 

FUNCTION OSUX ( CR tBKN ,E «NI ,X ) 

DIMENSION CR( iUOI tSKNI 1001 
SX ■ 0*0 
DO 65 l«l .NI 
BKNX m BKN ( I) **X 

SX ■ SX*(CRU »-EXP(-E/BKNXn«<EXPC-E/BKNX> >♦( AL06( BKN<I)))/BKNX 

65 CONTINUE 
OSUX « X^SX 
RETURN 

FUNCTION 0SUE(CRiBKN,E«N1 ,X ) 
dimension CR( 100 ) ,6KN( IOO) 

SE ■ 0*0 
00 66 1*1 ,NI 
BKNX a BKN ( I) ••X 

SE ■ SE + ( CR ( 1 1-EXP ( -E/BKNX »>•( EXPI -e/BKNX )) /BKNX 

66 CONTINUE 
□SUE a E*SE 
RETURN 

END 



PROGRAM CDCLEV 


This program computes free molecular drag and lift coefficients 
for flat surfaces at various angles of attacks and for a specified 
range of speed ratios. 



WKUN./TP cl>CLEV.UAHXXXXXXXXX,FANTASTIC-GK,5,I50 


InPLlCIT REAU*8( A-H.O-Z) 

01 ME MSI ON AN<i0(20)(CL(25|20)tC0(25*2n) 

REAOCi^tlOO) PJtALJtX 

Sw2«l .N1<4213B6237 

PI»3t l‘4159265i 

Su(PI«l .772H538b09 

S(.(2»l .9 1921356237 

E»2. 718281828959 

XK»X«PI/180. 

P«0. 32759 1 1 
A 1 «U, 259829592 
A2»-0. 289996736 
A3-1 .921913/91 
A9.-1 .953152027 
A5«l .061905929 
L0D«9.*3./P1 
00221 >1,25 
S*= I 

DO 22 J«l,10 
ANG-XR* ( J- 1 1 
ANGD ( J ) >ANG* 1 80. /P I 
SA«SIN( ANG) 

CA>COS ( ANG ) 

Z»S*5*.5*SA*SA 

1 -Z 1 


T>Z/3.75 
ES>E** ( -5*5 ) 

R> 1 .0/ ( 1 .0 + P*S ) 

tKFS«1.0-(Al*R>A2*R*R4A3*lR**3.)+A9«(R**9.)4>A5*(R**5.))*ES 
F*SQPI»(S*S*l.0-(0.25/(S*S)>l*ERFS*«S^«0.5/S))*ES 
COO-2. •F/(S*S*SQPl 1 
IF ( Z'3.75 MOi lOf I 1 

faU-(Z**.5»«EZ*( 1.0 + 3,5156229*T*T ♦3,0899929- (T«*9.) 

1 -1 .2067992* ( T--6. 1 ♦. 2659732* ( T6-8. ) ♦. 0360768* ( T** 1 0. | 
2*.0095813*(T**12. 1 ) 

Bl-(Z**1.51-EZ*(.5*.a7890599*T*T ♦.51998869*11**9.) 

2*» 15089939* ( T**6. I ♦•02658733*11**8,) ♦• 0030 1 532* 1 1** 1 0 • > 
3^. 0003291 1*(1**12. ) ) 

GO in 12 

0O-.39899226+.O1328592/1+. 0022531 9/ 

1- .00157565/1 1*1*1 )♦. 00 91628 1/1 1**9. 

2- . 02057706/1 1**5. )*.02635537/ll--6. 

3- . 01697633/(1**7. )♦. 00 J92377/ll-*8» 

B 1-. 39899228-. 03988029/1-. 00 3620 I 8/ 

l*.U016380l/(T*l*l)-.01031555/ll**9. 

2-. 028953 12/ (1**6.) *.01 787 65 9/ lT-*7. 

AL-8. -SA-SA/3. 

BL--9,-SA*SA*PI*SA 
C--Pl-.5*5A^(b.-SA-5A/6.) 


1*1 ) 


1*1 I 

♦ .02282967/ 1 1**5. ) 
-.00920059/1 T**8* ) 


Aj-ALJ 

Cl(I,j)-5Q2*SGPI*.5-CA*1(B0^B1 ) *A J* ( AL^BL*P J^C*PJ*P J ) 

1 ♦‘3 0*11. ♦13. -PJ)*AJ)/(S*S)^Bl*(l**ll ,/3.*PJ/3« )*AJ)/1S*S>> 
C<.(I,J)«(L00/(1.*L00))*CL(I,J) 

A0-- I . *9 ,*SA*SA/3. 
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Ty - 


aU«Pl-0.*5A*SA/3j 
CL>0».‘»,/2.-Pl ♦5.*SA*S*/i« 

0 0 • - . b ♦ P I / H . - S A • S A / 3 t 
K0»AD*BD*PJ + CC0*PJ*PJ-*t)0*PJ*P.J«Pv/ 

IF( J-l >20,21 ,20 

Clj(I,J)«SQ2*SA*((80*8n*(l. + Aj*K0>*B0*<l * ♦ A J • ( 3 • -2 , A;P,J> , 5*P J *P J ) ) 
l/l2t*S*S>-»Bl«(l. + ALJ*(l./3.*2.*PJ/3,-Pv>*PJ/<»t))/<2.*5*S)) 

2 + S(J2*(B0 + Bl)«(l.*ALJ*(-l# + l . 5*P J»PJ- , 5*P J*P J*P J ) )/(SA*S*S) 

CD( I ,J)-C0( I ,J)*SQPI 

COniJ)»«LOD/(i.-*-LOD))*CO(I,J>*COO/|l.+LOD) 

60 TO 22 

CO ( I , J » «SQP I • ( I . + A J* ( -I .♦ I .5«PJ*PJ-.5*PJ*PJ»PJ ) ) /S 
C0( I , J)»(LOO/( 1 .♦LOD) )*CD( I ,J>*CDO/( 1 .♦LOO) 

CONT 1 NUE 

WK I TE ( 6 , 200 ) P J , AU J 

WK I T E ( 6 • 20 1 ) ( ANGO ( J ) t J> 1 t 1 0 ) 

00 2 l«l ,2b 

WRI TE(6,202) I , (CD( 1 ,J) ,J«l , 10> 

WK 1 T£ ( 6 , 203 ) P J , AL J 
WKITE(6,201)(ANG0(J),J«1,I0) 

00 3 I « 1 , 25 - 

^RITE(6, 202)1, (CL(1,J),J>1, 10) 

60 TO 50 

FOWMAT ( 30 I 0.5 ) _ 

FORMAT! IHI ,25X,9HC6(S,ANG) i 0 X , 3hP J» , £ I 5 . 5 , R X , MH AL J» , E I 5 • 5/ ) 

FORMAT ( HX , 7HS. . ANG« , loe 1 2 .5/ ) 

FORMAT ( 3X , 1 2 ,6x , 1 OE 1 2.5 ) 

F0RMAT(1HI,25X,9HCL(S,ANG) , 8 X , 3hP J« , E 1 S • 5 , R X , RH AL J- , E 1 5 • 5/ ) 

END 
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PROGRAM CFEVAL 


This program computes values of force coefficient for a 
specified speed ratio as a function of gas surface Interaction para- 
meters and angles of attack. 
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6RUN ./TP CFEVAL .UAHXXXXXXXXX .FANTAST IC-GK «B » IBO 


implicit REAL*ei A-H.O-Z » 

dimension AI50.2) .6(B0|2) .HIBOtS) iQ(50i 2) tP(50t2) 

DIMENSION E(50«2) 

DIMENSION GAM(20) ,ALJ( IS) .PJ( IS) .CF|20,20,20) 

DIMENSION AS(BU,2) |GS(S0.2) .HS(SO,2) ,QS(50,2) .PSI50.2) 

DIMENSION FF(50,2) 

READIS, lOOISRAl 
SRA2-SRA 1 *SRA 1 
AR-l .5780 
A 

A2«.06260A01 220 
A3«.0‘47573835H6 
AH-.Ol 736506M51 
Bl>. 24998368310 
B2«.0920U180037 
B3» . 04069697S26 
BH- .00B26449639 
AAO- 1 . 38629436 1 1 2 
AA 1 -0.09666344259 
AA2-0.0359U092383 
AA3-0. 037425637 1 3 
AA4-0. 01451 196212 
bOU-O . 5 

681-0, 12498593597 
BB2-0. 06880248576 
BB3-Q. 03328355346 
BB4-0.0044 1 7870 1 2 
PI-3. 14159265359 
U051 1-1,19 

CiAM( 1 )-5,0*( I-l ) 

5AMR-5.0*U-n»Pl/l80.0 
5 1 NG-S I N ( GAMR ) 

COSA-COS ( GAMR ) 

SlNA-SING 

FM 1 - 1 , 0-S 1 NQ*S I NG 

FM2-FM I *FM 1 

FM3-FM1-FM2 

FM4-FM1-FM3 

J- 1 

IF(FMl)21,31,21 

t<I.j)-l.+Al*FMI+A2*FM2*A3*FM3 + A4*FM4'«-<Bl*FMl+B2*FM2'«>B3*FM3 + B4*FM4 
1 ) *ALn& ( 1 ,/FMl ) 

FF11,J)-aa0 + AAI *FM1 ♦AA2*FM2*AA3*FM34.AA4«FM44.1BB0^BbI*FMI*BB2*FM2 
^♦BB3•FM3■^B84#FM4)*U06( l./FMl ) 

GO TO 41 
£ ( I . J ) - I . 

FF « I , J ) - 1 000. 

A( I ,J)«?.+4.*AK*E< I ,J)/PI 
ASl I ,J)-2.^5,*AK*FFn iJ)/PI 

G1I,J)«AK*4,»E(1,J)/(3.*P1 ) -fAR- 1 6 .*cOSA»COSA»FF ( I ,J)/(9,*PI ) 

GS( 1 , J)--AR*FF« I ,J)/PI 

H(I,J)-4./3.+Ah-(EtI,J)*(-88./(9.*Pl)*4,-16,*C0SA*C0SA/(9.*Pl)) 
l+PI*SINA*SINA*.5-8.*FFll,J) -COSA-COSA/ ( 3. »P I ) ) 
mS(1,J)«4./3.+AR*2.*E(I,J)/(3.*p1)4aR*FF(I,J)*(-1.+2./(3.*PI)) 
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J+AR*(-Pl*SINA*SINA/<<»-PI»3.*(SlNA«»i*.)/l6.) 

«( I iJ)»l #/6.+AR«E« I iJ)*( 16./(3.*PI I -*♦ . *8. •COSA*C.QSA/,( 3. *P I ) ) 
I+Ah*(6.*E(I ,J)/PI-.25*P1«SINA*S!NA + 8**C0SA*C0SA»FF( I , J ) / ( 9 . *P I ) ) 
QS( I ,jl»l./6.-AR*E( I |J)*2«/(3.*PI )+aR*FF( I 1« + 16#/(3«*PI ) ) 

l-AR>1.5*FF(I|J)/PI^AR«PI«(SlNA»SlNA«(l*/8*-^3t*SlNA»SlNA/32*)) 
P(I,J)--1./6.+AR*E(I,J)*(-8./(9**PI)*1.-8.*C0SA*C0SA/(9»*PI)) 
1-AR*2.»E( 1 , JI/PI 

Pb{I,J)»-l./6.^AR*EUiJ)/(3.*Pn-AR*FF(IiJ)*(.25 + 2./PI) 
l+.b*FF( I ,J)*AR/P1 
DO bl J«1 , 1 1 
UU bl 1 , 1 1 
ALJ< J)».2*( J-l ) 

PJ(K)»#2*(K-1 ) 

PJ1»PJ(K ) 

P J2»P J 1 «PJ 1 
P j3«PJ 1 *P J2 

CF(K,J,n«A(I,l)*AS(I,l)/SRA2+ALJ(J)*(Q(I|l)*PJl*H(I,l)>PJ2*0(Iti) 

l+Pj3*P(I,l))+ALJ(J)*(6Sn,l)+PJl*HS(I,n+PJ2*QS(I,n*Pj3*PS(I,l») 

2/SKA2 

CF 1K,J,I )*CF(k;iJi 1 )/J l.+AR) 
bl CONTINUE 

00 lb J« 1 . 1 1 
WRITE|6,2I 1 ) JfAUJ(J) 

21i FORMAT! 1X,**HALJ( I I2,**H) ■ F8.2) 

WRITE! 6,9797 ) 

WH 1 TE ( 6,9997 ) 

WR 1 TE ! 6 , 9899 ) 

00 1 1 IS K»1 , 1 I 

WRITE ! 6,21 1 1 > ICF!K , J, I ) , 1-1 ,7 ) 

WRI TE I 6 ,2 I 1 I ) !CF ! K I J, I ) , I-8| IH) 

WRITE!6,21in!CF!KiJ,I)il-l5,19) 

111b WR 1 TE ! 6 • 1 62 ) 
lb WK I TE I 6 , 1 63 ) 

GC TO 99 

1 62 format {//) 

163 FORMAT!////) 

2111 F0RMAT!IXE15.9,2XE15»9|2XE15.9»2XE15#9»2XE15,9.2XE15,9,2XEI5.9) 
9797 F0RMAT!7xtlH0il6XilH5,16Xf2H10tlbX,2Hl5,lbX.2H20tlbX,2H2b, 1BX.2H30 
1 ) 

9997 F0RMAT!7Xi2H3bfl5Xt2H9U,15X,2H9b,lbXi2HB0, IbX ,2Hbb t 1 bX , 2H60 , 1 bX • 2H 
16b ) 

9899 format! 7Xt2H70t lbXt2H7b,lSX,2H8Q, l5xi2H8B, 15X|2H90) 

100 F0RMAT!D10.b) 

ENO 


132 



PROGRAM RUFSPH 


This program was used to compute drag coefficients for a 
number of non-spherlcally shaped objects. The copy shown here was 
used to conq>ute the drag coefficient of ellipsoids of various eccentric- 
ities and gas surface Interaction parameters. 
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WKUN*/TP RUFSPH,UaHXXXXXXXXX,ENVIR-DYN-GK, 5, IBa 


implicit REAL*8(a-HiO-Z» 

Dimension iThbizi > tTHp(2n »txbi2i i ,txfI2i) 

DIMENSION DI (21 ) »0R(2| *21 ) ,X(2l I ,W(21 ) *P<J(21 I »COF(2l »2l I .6 AM (2| > . 
ITHI (21 ) fWW(21 ) 

DIMENSION THX(NO» 
dimension qA( 10) iGA2( 10) 

DIMENSION GR(20) 

DIMENSION F I ( 21 ) ,SS( 21 ) 

DIMENSION RR( 21 ) 

X( 1 )>0*07AS2AB21 I33B0 
X(2)«0.2277858B1 IHUB 
X(3)-0.37370608871SH2 
X ( ■0.51086 700195083 
X ( 5 ) ■0.63605368072652 
X ( 6 ) ■0.7H6331906M601B 
X( 7 )^0«8391 1697182222 
X(8)«0.91223H<42825133 
X(9)^0. 96397192727791 
X( (□(■0.993128B99185I 
W( 1 )■0• 1527533871307 
W( 2 ) ■□. 1 H9 1729869726 
W( 3)^0. 1920961093183 
W( >M^0. 1316886389992 
W(5)«0. 1 18199B31961B 
W(6)»0* 1019301 198172 
W(7)^0.0B32767915767 
W(8)^Q. 0626720983391 
W( 9 ) ■O. 09060 19 298009 
W( 10)«0t017619007|32 
P^O. 327591 1 
A !■□. 259829592 
A2*-0. 289996736 
A3^l .H21913791 
A9--1 .953152027 
A5^1. 061905929 
PI^3. 19159265358976 
DO 2 !■! , 10 
l-I 

GAM( I )^PI*.5»( I .-X( J) ) 

2 WW( I ) J ) 

0031^1 ,10 
J^l^lO 

QAM( J)^PI*.5«( 1 .♦XI I ) ) 

3 WW( J ) »W( I ) 

DO 90 M^ 1 ,S 
ECC ■ 0.0001*M 

GA2(M)^.5*( 10.««(M-1 ) )*PI/( 360 ••3600. ) 

GA ( M ) mECC 
AS^l .0 

BS»( I .-ECC*ECC)«*.5 
D09M1 ,21 
009J^i ,21 
9 DR( I , J)«0.0 
DO 29 !■! , 10 
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S(,A»SIN(6/Ii+1( I >"1 
C(,A>C0S(GAM( I ) ) 
C2>CGA*CGA 


S2"SGA*SGA 

RR( I )>(C2/(BS*BS>-«-S2/( AS*AS) 

DR0TH*( (RR< I ) l»*3* )«(BS**(*2* )•AS•*(•2• ) )«C6A«S6A 
GR( I )«ATAN(0R0TH) 

29 continue 

DO 27 I>1 1 ,20 
GR( I )«0,0 
RR( 1 )>t .0 
27 CONTINUE 
00SI*1 • 10 
SS( I )b10,«»(.1 ) 

T«J ./( 1 .♦P«SS( I ) ) 

S>SS ( I ) 

EKF5«2t*(S-S«S*S/3.+(S**5.»/l0.*(S»*7.)/M2» 
l*(S**9. )/216. I/(PI**.5» 

01(1 )m0.0 
006 Ja 1 , 20 

Y-ABS ( S»COS ( GAM ( J ) ) ) 

T«1 ./( 1 .♦P*Y) 

ERFYa2.*(Y-Y*Y*Y/3.*(Y*a5#l/IO.-(Y**7t)/M2. 

!♦( Y**9. )/2l6. )/(PI**.5) 

CGAaC0S(GAM( J) ) 

1F(CGA)31 ,33,33 
31 Ya-Y 

ERF Ya-ERFY 
33 PQa,5*,5*ERFY 

QaS*Ya(l.+ERFY)+5aEXP(-YaY)/|PI*a,5) 
Fl(J)a(Q«q«2.»PQaQ«CGA>PQ«PQ)a*iS 
THI ( J)aPia,5.AC0S( (q*CGA*PQJ/FI ( J> ) 

THX( J)aPla,5-THI ( J) 

6 CONTINUE 

00 70 Jal ,20 
ANB-GAM ( J ) -GR ( J ) 

ANF aGAM { J J -GR ( J ) 

YBaABS ( S*COS ( ANB ) ) 

YF»ABS( S*C0S( ANF I ) 

TBa 1 ,/ ( 1 .♦pays ) 

TFai ,/( I ) 


ERFYBa2.*(YB-YB*YB*YB/3.a(YB*»5,l/lo.-(YBa*7#)/H2.+(YBaa9.)/216*)/ 
I ( PI ••.5 > 

ERFYFa2.*(YF-YF*YF»YF/3. + (YF**5,)/Io»*'<YFaa7.)/R2«-f(YF«*9,)/2l6.)/ 
I (Pl**.5» 

CANBaCOS ( ANB ) 

CANFaCOS ( ANF ) 

IF (CANB >51 ,B3,S3 
51 YBa-YB 

ERF YBa-ERFYB 
53 POBa.Ba.saERFYB 


6 1 


QBaS*YB*(l«*ERFYB>*SaEXP(-YB*YB»/(Pl**.5) 
FB( J)a( UB*qB-f2««PgB«QB«CANB<fPQB«PQB)***S 
THB(J)aPIa,5.AC0S((QB«CANB«PQBt/FB(JI> 

IF (CANF )61 ,63,63 
YFa-YF 

ERF YFa-ERFYF 
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63 PQF»#5*#5»ERFYF 

QF«S*YF*( 1#*ERFYF j*S*EXP(-YF*YFf/(Pl**#5) 

FF ( J »■! QF*QF*2.*PQF*QF*CANF»PQF*PQF »»**5 
THFU»*PI*«5-AC0S( (iaF*CANF*PQF)/fFljn , 

TXB< J»«PI«*5-THB( J) 

TXF( J).PM*5-THF( J| 

70 continue 

OO 8S L> 1 i20 

OiB«fB(U*COS(gAH(L i-Pl**5«THB(L)*GRU.) > 

|•RR(L)•RR(L)•WW(L)*SIN(6AM(U )/C0S(QR(U ) 

OKI )>OIB*OI ( I ) 

85 CONTINUE 

on I )«oi ( n«Pi*tS 

007K>I (21 
PJ(K)«. 1»(K-I I 
00 8 L> 1(20 

THRF»PJ( X J^PK.B-t’ncPJIKn^THFtU 

0RF“FF(L»«C0S( THRF*PI*t5-SAMU.KGR(L) )*WWID*SIN 
1 (6AM(L) )*RR(L) *88(1 ) /cos (GRID 1 
0R( I ,K )«0R( 1 (K l-ORF 

8 CONTINUE 

DR( I |K )«PI««5*0R( I iK ) 

coF ( I (K ) « 0 R( I (K ) /on n 

7 CONTINUE 
5 CONTINUE 
009I-1 ,10 

MRiTE(6, ion SSI 11 fOt ( n 

WRITE(6,1031 GAIN) 

OOlOKsl ,21 

WRITE|6,t02) PJIK),C0FII,K),DR11,K) 

10 CONTINUE 

9 CONTINUE 
90 CONTINUE 

WRITEI6, 105) 

DO 28 1-1,20 
AI-I 

GR( 1 )>GR( I )*180./P1 

GANI I )-GANI I )«180*/PI 

WRITE16, 10H)AI ,GAMl I ) , RR I K , GR I I ) 

28 CONTINUE 

105 F0RMATI//9X,lHl,6x,6HGAN|n ,6X ,5HRR ( I ) ,6X ,5HGRI I ) / ) 

109 F0RMATI2X,F10.5,2X,E10*5,2X,E104S,2x,E10*S) 

103 F0RNAtI9X,3HGA>, £10,5/1 

101 FORMATI ////9X,2HS-,E10,5,2X,6H0I I S ) ■ , E 1 3 • 5/ /8x , 2HP J , 6X , 7HK ( P J , S ) | 
|5X,8H0RIPJ,S)/) 

102 F0RMATI9X,F10*5,2X,E10.5,9X,E10.5) 

200 FORMATI / / / /9X , 2HS- , E 1 0 • 5 , 2X , 6HD I IS)«»E10,5,//8X,1HK, lOX ,6HGAMI K ) « 
16X,6HTHI(K ) ,6X ,5HFI IK ) / ) 

201 F0RMATI9X,F10*5,2x,El0,5,2X,EI0*5,2XtEl0*S) 

END 
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PROGRAM CLELAM 


This program was employed to calculate orbit perturbations 
due to lifting satellite shapes. 
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n n n r\ rt r\ rt 


9ttUN,/TKC CLELANiUAHKXXXXXXXX tORBiTtlOtiOO 


C calculation of SATELLITI 0 R 8 *T change caused by lift and drag 
c PENTURbaYIONS..... NUMEBICAL integration using a NODIFICD RUNga* 

c kutta technioue 

c 

implicit double precision (a*H« 0 «*Z) 

dimension YI7),«(7)tQ(7).PHlNT(IO) i PH 1 IM 1 01 i A A ( 5 ) » Y Y ( 7 *40 ) » OeN ( AO 
DtHTIAO) 

double precision mu 
real ESTEPiOIFFEQiRtO 
external OIFPEO 

common rAiMUtP,PF*ACT«NtR«ACN»AtK«VMEAN 
INPUT OF INITIAL VALUES 

L ■ number of ECCENTRiCf T ies times number of perigee heights 
M ■ number of Different satellite attitqdes 
N ■ number of area/mass Ratios 

NO - upper limit on The number of orbits for a single launch 
READ (S| 3 ) N,M«L(N 0 
3 FORMAT (MIS) 

READ (StM ) ( AA ( I ) ( IP I fN ) 

M FORMAT(EIStT) 

READ ( 5 , IM) (PHINT(K) ,PH 1 R()(| , K ■ 1 ,M) 

IM format ( 2 F 1 S.I 0 ) 

REA 0 ( 5 . I 0 )VE.RE,MU 
10 F 0 RMAT( 3 EIS«Y) 

READ (S,I 3 ) GA*PJ 
13 FORMAT ( 2 FI 0 * 6 ) 

READ (S,| 5 ) ( (YY(MM>NN| ,HN • 1 > 7 ) .OEN f NN ) • HT ( NN ) t NN « IfU 

IS FORMAT ( 5 EI 5 « 9 / 2 E 15 « 9 / 2 tlSt 4 ) 

Pi ■ 3. IM15926S359 
PSTEP-.2 

STEPEXP , 3010299954 

DO 1010 I ■ I iN 
AM ■ AA ( I ) 

00 1020 K ■ I »M 
PI ■ (Pl/I 60 t ) • PHlNTlK) 

P 2 ■ (Pl/| 80 t) • PHI* <RI 
*R|TE ( 4 fl 2 ) NfMfLiNO 

12 FORMAT(|H|t *FLaT PLATE IN aN ELLIPTIC ORBIT ABOUT THE EARTH*, 1 SX, 
IMIA,///) 

*RITE(At> 4 ) AM,PHIW(K ) «PHINT(K) , 6 A,PJtVE»RE» MU, ESTEP 

16 format nx, ‘Parameters' ;///» * satellite area/mass ■ » » isx»ft.3,9x 
I » 'M2/K6' »//i • ATTITUDE ANGLE PHIR ■ • » I SX *F7 • 3 #9X , * DEG * » // t • ATTiT 
2U0E ANGLE PHINT a t , I M X ; F? • 3 « 9 X , « DEG » t /> » * GA (GAS SURFACE PARAMET 
3ER) ■ » ,8X iFfl*M ,20X , ♦ CGaMMA) X SORT I I -ALPHA ) REF* I'*//* 

M» PJ (GAS SURFACE PARAMETER) • * i 8X i FB . M » / / , • EARTH ANGULAR VELOCI 
sty ■ *illX,ElM*8* IX»»RAO/SEC',//,* EARTH RADIUS • 'iZSKiEtOtG 
6i3X.*Mt;//t* earth gravitational CONSTANT m * »7X «E U *7*2X1 'HS/SC 
7C2',//,' ECCENTRIC ANOMILY STEP SUE • 1 » 7X , F 1 1 .7 , 5X * * RaO • * /• 1 H I ) 
AT«C0S(P2)*C0S(P1 ) 

AN-ATpTaNIPI ) 
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II nil 


AK»S|N(P2) 


H 


maim orbital CLEMENTS ' 

Ti-eeccNTmc anomaly 
T2«SEMIMAJ0II axis 
Y 3»CeCCNTRIClTy 
TRaORBITAL INCLINATION" 

YS-aNCUMCNT or PCNiaEC IOITM NESPECT to A5CEND1N6 Node) 
Y6-L0NSITU0E Of - ASCCNDttICi NODE f«iTN NESPECT TO VERNAL CQUJNOX) 
T7-TIME or RERISEE PAlSBiE ■ - 

TIMCmTIME INTO ORBITIS) 


DO 1000 I 1 ■ I ,L 
00 20 MM ■ I ,7 
20 Y(MM) ■ YYIHMt I I ) 

RHOO ■ DENI I n 
HH - HTI I 1 I 

HPI ■ Y(2I«(I •«Y(3) )-RE 

VP*SQRT< IMU/YI2) )•( I f»Y(3) )/f t t»Yf3) ) ) 

DYNAPP« RH00»VP«*2 

0LDY-0«0 

0LDT»0«0 

norbit^i 

ITER-0 

J • I 

NULL *0 

CALCULATION Or VARIaOLES 

2<« NEXTbO 
NOIP • 0 
25 1TER«ITCR«^I 

R ■ Y(2)*( I •«Y(3)*C0SIYn ) ) ) 

HmR-RC 

RP-YI2)*! I ••YI3) ) 

HP ■ RP • RE 

ir fHPtLT«0*0*0R*Y(3) (LYtOtO) 60 TO VM 
RH0>RN00 «EXPI-(wNPI )/HH) 

rr - SORT! I 1 .-Y I 3) ••21*1 1 •*YI3)*C0S|Yt 1 I ) >/I 1 • *Y I 3 ) *C0S f Y | | ) ) )) 
P ■ YI2)*I I ••Yf 3)«*2) 

TA - ACOSI tCOSI YC 11 )-Y(3) )/f I ••Yf3)«C0S(Y| 1 ) ) ) ) 

TaSIN • tSINIYf l))*SQRTM«*Y(3)**2))/tl»-Y(3)*C0S(Yll))) 

IF (TASINI H0f5Q,50 
NO TA»(2.»PI I-TA 
50 U«TA*Y(5) 

AR>| 1 t/Fr)*IAT*Y(3|«SlN(TA)»AN«(tt*Y(3)*COS|TA) ) ) 

AS«( 1 «/rr )«IAT«( 1 •4Y(31«C0SITA) )»AN«Y(3)*S1N(TA) > 

AaSQRTI AR**24AS«»24AK*«2) 

VR- SQRT(MU/Y(2) )*tY(3)fSiNfYMU)/(l«-Yl3)»C0S(Y( I ) ) ) 
VS>SQRTIHU*( l•4Y(3|•C0S(TA) }/R)*(VC*R*COS(Y(R) ) ) 
VK»VC*R«COS(U)*cOS(Yi<l)T 
VmSQRTI VR**24VS**24VKP«2) 
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VMCAN»S«RT(HU/r (2)**3» 

DOT«(AR*VR^AS*VS«AK«VK)/( A«V) 

IP (ABS(OOT)'fLT.ttO) 60 To 51 
OOT«1 . 

51 AN6ATK<iAC0S(D0f) 

TAS*AN6ATK 
IF (TA5) 55i5*«,55 

SH COTAS ■ 0*0 
60 TO 5* 

55 COTAS • COTAM (TAS) 

56 TAJ ■ IP|/2.*TAS) 

CDV«2,*SINC TAJI w2,*QA*SlNlTAJ»*C0Sf (Pl/2» I *P J* ( 2 J > *T A J > 

COR«(VR/V)*COV 

COS*l VS/VI*C0V 

COK«( VK/V J*COV 

CLZ»*2«»6A*S1NITA J)«S1N1P|/2»«PJ*(2*"PJ>*TAJ1 
CLR«( AR**I VR/V)*SINCTAJI I *CLZ/S ! N | TaS ) 
CLS»(ASP<V5/V)#5lN(TAJ)»*CI-2/SlNtTAS) 

CLK«(AK«IVK/V)*SIN(TAJI » *CLZ/S I N I TAS I 
CDV ■ COV*l •R867T2I 17 
CtZ ■ CDV«C JI/SINIR2I » 

AHM ■ AR/2*00 
RL ■ CLZ*AMM 
BO ■ CDV*AMM 

ACR«(RhO*V**2I*( CBL •AR/( VMEAN*SINITASM >•( I ••YC3)*C0SCYI I ) ) >»Y(2 
I )«Y(3>»SIN(Y( I ) )«(BL • COTaS *80 >/V) 

ACS«IRH0*V**2)»I (BL *68/ I VRE AN«S I N ( TaSIM • «l •«»Y 1 3 > *C0S I Y C I )» ) • 

I (SQRT(i;-Y(3I**2>»(VE*C0S(YlRI>/ VMEAN)*f|»-Y(3)«C0SfYnni ♦•21* 

I (Y(2)/V)*(BL • COTAS *B0 M 

ACK»<RH0*V)*I (VE*Y(2»/VMEANI« SINIYIDI *€05 1 U I • C (I . -Y ( 3 1 • COSCYC 
lini**2) •I-BL vCOTAS •bo t*fBt •V/VHE|lNl*AK»llt»Y(3)*C0S<Y(|)m 
ACN«( I #/FF»*I ( I .♦Y(3»*C0S(TAl I • ACR* Y ( 3 1 *51 N ( TA 1 • ACS) 

ACT«( I ./FFI*f Yf 3I*S1NITA1«ACR** 1 1 *Y C 3 1 *£08.1 TA ) 1 *ACS 1 
TlME*CVn )-Y(3»*SINIT( 1 tl 1/VMEAN ♦Yf 7) 

RA"Y ( 2) • ( I • ♦Y ( 3 ) ) 

OYNA ■ RH0»V**2 

If I(0 Y»IaP0YNAPP|.UT*I*00*5) 60 TO 58 

ESTEP ■ IPSTeP*fOYNA/OVi|App1**l*STtPEX> »**0I7R532»25I 
GO TO 5? 

58 ESTEPa* I 396263R0I 

59 ORBITpNORBIT 

IF (NOIP«EQ«l tOR«NEXY,EO«|) 60 TO 81 
IF ( NORBIT.EO. 1 • AND* ITERtEQ* 1 I 60 TO 60 
60 TO 70 
60 OLDRA ■ RA 
OLOT ■ TIME 

tai-ta 

Y2«Y<2) 

Y3*Y(3) 

YR-YIRI 

Y5«Y(5) 

Y6«Y(6) 

HP1»hP 

0L0HP»HPI 

0LDY3-Y3 

V1«V 
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I 


ACTHiACT 

write fi.lS) (YfMN) tHN"l t7) •TA.HRtKHOO 
18 FORMaTIJBX, ••••••••*• lillTlAL MAIN ORBITAL CLfcHENTS ••••••••• t ,//^ 

1 ilX,tECCENTR1C SEMIMAJOR ECCENTRICITY ORBITAL ARGUMENT 

20f longitude Of*. ax. ‘T ime of*»8X.»true«»sx, •perigee Rerigee*./. 

3* ANOMALY*t4X|«AXlS*.2B)(,«|NCLINATlON PERIGEE ASCENDING NOd 
he perigee Passage anomalt height density • •/•rx, *yi * , |0x, 
5*Y2* . I IX) *Y3* *12Xt«Y<«* .l2X.*YS*.12Xt»Y4» .I5X»*Y7* .|2X»*TA*tBX»*HP* 
A.7X ( *RHOO' .//.Ff •A. e 1R«7« p9 . R .F 1 4 • 7 . 2F1 R • 7 .p] 7 .7 »p 1 R .S .XeI 0 • R . / • / 
7.B2X, MM-K«S-RaOI ANSI • ;//// I 
70 IF < Y( I I-C2.RPI I I 70.80.80 

80 NORBiTiNbRBitil 

81 oelra«oldra-ra 
oelt.time-olot 
0CLTA*2.*PI*ITA-Ta1) 

DELA-Y2-YI2I 

DELACT-ACTl-ACT 

0ELV«V|-V 

OELHP*HP*HPI 

DELC«Y<3»*Y3 

0ELI«Y(R»^YR 

0ELR*T(5)-Y5 ]] 

0ELRR»Y(4)-Y4 , 

ADECAY ■ OELRA/TIME 
THETA*I80.»TA/Pi 
EC*I80.*Y( I )/PI 

IP I NDlP.EQ.l.OR.NCXT.EQ.li GO TO 420 

OLORA - ra 

olot ■time 
tai-ta, 

Y2*Y(2» , 

Y3«Y(3) 

YR»Y (R ) 

Y5»YIBI 
Y4«Y ( A » 

HPl-HP 

oldhp-hpi 

0L0Y3-Y3 

V1«V 

acti-act 

IF (NULL.EQ* I I GO TO 82 
Y( I )«Y ( I )-2.*Pl 

82 DYNAP«(RH0*V**2)*.0000l 
0YNAPP»(RH0*V**2) 

OUTPUT 

1F(HP.6T.0.0*AND.Y(3) .GT.0.0 I GO TO 410 
write (&*700) 

700 format i//,iox,t Satellite has crashed on this orbit or The orbit 
I HAS been reduced TO C I RCULAR • . /// I 
410 WRITE (4.8001 

IF (NULLtEO. 1 I 60 TO 420 
write (4.8011 

420 WRITE (4.810) NORB I T . THETA .EE »V ( 2 ) . V ( 3 ) .Y 1 R ) 

WRITE (4.82QI 

WRITE (4.830) Y ( 5 ) . Y ( 4 ) . Y < 7 | . T | ME 
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xj \j %j nj %j %j 


,111 II nil I 


III IMIHIIIIIII II IHIIIHHMHIl I ■ I I 


write (6tSR0) 

WRITE U.tSO) VRiVSfVKty iACN|ACT,ACK 
WRITE (AtSAO) 

WRITE <*,S70) R,H,HWiW»RHO.TAS 
WRITE (AftSO) ADCCAY 

WRITE IAi«aS) DEtHP»OELttOELl iDELWiOELWW 

WRITE! A,SWO)OELTA,DELA»OEURAtDei-T,OELV»ClELACT,ITER 

WRITE !At7»S) CL2 .CU^ »Ct* ,ClR iCOV , CDR . COS . CDK » A » AR » AS f A« 

IF (NULLtEQal) fiO TO 1000 
IF (NOIPtEO* i ) 60 TO 2R 

IF (NEXT«EQ* U 60 to Z** 

795 format (37X,*R*,l6X,*S«iISX,*K*,/ ,* LIFT COEFP I C I ENT » , / ,9X ,90 1 9 1 
no./ ,» 0RA6 coefficient* ,/,9X, 9019. 10,/ ,» ATTITUDE VECTOR*. /, 
|9X,90I9,|0»//////I 

fiOO format ( IX, *NORBTT*;i 1X,*TA*,I7X,*YI*,17X,*Y2*,I7X,*T3*»17X,*Y9*,I 
IOX,*TA,Y| IN decrees*) 

801 format ! I X , *NEAR FeRICEE*) 

810 format ( I5,9X,SDI9.|2,/T 

820 format (I8X,*Y§i,17X.*T|*,17X,*Y7*,|AX»*TIME*) 

830 format (9X, 9019, 12,/) 

890 FORMAT ( 8X * * VR • , 1 A X , * VS * » I AX . * VK * . 1 A X , * V * . I 6X , * ACn * , I SX , * ACT * , I 5 X , 
I *ACK*) 

850 format <7018.1 I ,/) 

8 AO format (9X,*R*,l9X,*H*.18X.*HF*,t9X,*RA*.18X,*RH0»,l8X,*TAS*) 

870 FORMAT (AD20.13,//) 

880 format <AX,* THE AR06CE DECAY RATE • *. 1020.9, » MeTeRS/SECONO • , / / ) 
885 format l7X,*DELHR*.l3Xi«DCLE*.l9X,*8ELl*.lSX,fDELW*,13X,t0ELWW*./, 
15EI8.I2,/ ) 

8 90 FORMAT !7X.*0ELTa* , I 3X , * DEL A * . 1 9 X , * DELR A * . 1 9 X , * OELT * , I 9X , * OELV ♦ , I 3X 
I f »OELACT*,10X,*|TERATlORf* ,/,AEl8tl2.8X,I5./» 

ITER-0 
J • 0 

■ 0BBan0BBaR0DnB«iiaaa0)iRaBaiiBRaoRRBaaBBBBBaiinBBBBBnHBaBBaaaB*aoo0iB 

90 CALL 6ILL!0IFFE8.Y,.1 ,ESTeF,W,Q,7 ) 

bbbbbbbbbbbbbbbbBbbbbbbbbbbbBbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbb 

IF <HF.lt. 0.0. 0R«Y(3) .LE.0.0) 60 TO 999 

IF <N0RB1T»N0) 910,1000.1000 

910 IF <HF*6T. < *8«0LDHF) *AN0.Y<3) .6T. < *8*0LDY3) ) 60 TO 915 
OLDHFbHF 
0LDY3«Y(3) 

NOIF-I 

WRITE <A.800) 

60 TO 25 

915 1F<0YNA46T,0YNAF.0R. J.EQ.I ) 60 TO 29 

920 Y< I )«2,»FI*Y< 1 ) 

J ■ I 
next - 1 
WRITE <A,800) 

WRITE <A,802) 

802 format <1X,*NEAR AF06EE*) 

80 TO 25 
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rtooo nnnnrkCtnnn ® r»r»or>r»r» 


ff? NULL ■ I 

GO TO il 
1000 CONTINUE 
1020 CONTINUE 
1010 CONTINUE 
END 

OFOR.IS OIFFEQtOlFFEQ 

l^EAL function DIFFEO(Ytl) 
dimension Till 

double precision r tTAtHUtP.rr»ACTtU|RtACN»ACKtVMEAN 
COMMON TA|MU,P,FF ,act,u,r,acn.ack,vmean 

aoRnHnii0MaaaaBaaBkBnHiiBaaaBBa ■■■■■■»■ ■■■■■BBaRHaBaoBonBaoBliBBaHa.a 

differential equations for OMBITaL PaNAMCTENS,Y2^Y7, RITH RESPEcT 
TO THE independent PARAHETERt Yl 

aaBaBBRBaaRaaBBaBaBaaaaBBBBaBBBBBBBBaBaBBBaaaBBBBBBBBBBBaaBBBaBBB 
GO TO (3IOi320«330t3HOiSSO,3AO) (1 
310 01FFE0"<2»*<Y(2J**2J/SQRT(MU*P) >*FF*ACT 

60 TO ROO 

320 OIFFEQ-ISQRTI P/MU ) /FF ) • ( 2 . • I COS ( Ta ) ♦ Y ( 3 » ) • ACT*R* ( t . -Y ( 3 ) **2 ) *S I N ( 
ITa)«ACN/P) 

60 TO «I00 

330 01FFEQ»R*C0SIU)*ACK/SqRT(MU»P) 

GO TO HOO 

3M0 D1FFE0«(S0RT1P/mU)/Y13) )• I I 2 • *S 1 N I T A ) • ACT» I R/P ) • C 2 . *Y 1 3 ) *COS I T A ) ♦ 
I(Y(3)**2)*C0S(Ta) )»aCN»/FF*R*Y(3)»SIN(U)*ACT/(P«TAN( Y(i») » ) > 

IF (YIMH *IOO,3f9,NOO 

3S0 0IFFE0»(R«5IN1U)*ACK)/<$0RT(P»MU)«SIN(Y(*»M) 

IF (YIR)) MOO»3»T,««00 

3*0 OIFFEO-1 Y(21/MU)*(2.*(P/Y(3)*R»Y(3) ) *S I N ( T A ) /FF-3 » *SORT ( HU/P > •FF* < 
IYM)-Y(3)*SlN(Ynn )/VMEANl«ACT*nR*P*C0SITAn/(MU*Y<3l*FFn»ACN 
GO TO HOO 
3»Y DIFFEQ«OtO 
Roo Return 

END 

FORtlS SILL, GILL 

subroutine gill (DY,Y,Z,H,W,0,NI 
DIMENSION Y(N),W(N),QIN),A<<*),Ct*n,BIR) 

double precision Y 

DATA(All),C(n,B(I),lBl ,N)/2«*5,2, ,2*.2f28T3283, 1 • , 2* 1 t 707 J 0*7 1 , 
lit, • l*6&*****, •S,2«/ 

This routine is a modified runga-kutta numerical integration 
technique 


ox» IS The interval size. 

R- IS THE ARRAY USED TO STORE THE 
VALUE OF Y«(X)> WtllaFOIX)*! 


DX«H 
W( l)Bl • 


FOR THE FIRST INTERVAL THE Q*S ARE. SET TO ZERO* FOR 

subsequent intervals the previously computed Q*S 
ARE used* 
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rt r» r» 


C 

DO 5 J>1 tN 

I e<ji"0* 

to 00 20 J«l t*> 

00 IS K«2«N 
IS «(K)*0V(V»K«l ) 

00 20. K«1 tN 

Y(K)«VIK)«0X*A( J)*OtKn 
10 0(K)«Q(K)«3**A( Kl*8( J)«Q(K) }•CIJ)•0(K) 


test if value of inoependemt Variable 
HAS BEEN reached* 


1F( YM >♦DX•GT•Z) DX*Z«YM) 
IF(OXfLT*ABS( Y( 1 n«2*E<’8) 60 TO 2S 
IF(Y(|)-Z| l0t2S«25 
IS RETURN 
END 


*U.S. GOVERNMENT PRINTING OFFICE; 1976 - 735-004/13 
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